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Abstract 
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Nonequilibrium Green's function methods allow for an intrinsically consistent description of 
the evolution of quantal many-body body systems, with inclusion of different types of correla- 
tions. In this paper, we focus on the practical developments needed to build a Green's function 
Q methodology for nuclear reactions. We start out by considering symmetric collisions of slabs in one 

^ dimension within the mean-field approximation. We concentrate on two issues of importance for 

I actual reaction simulations. First, the preparation of the initial state within the same methodology 

If^ as for the reaction dynamics is demonstrated by an adiabatic switching on of the mean-field inter- 

T 1 

p\| action, which leads to the mean-field ground state. Second, the importance of the Green's function 

o 

matrix-elements far away from the spatial diagonal is analyzed by a suitable suppression process 
that does not significantly affect the evolution of the elements close to the diagonal. The relative 
lack of importance of the far-away elements is tied to system expansion. We also examine the 
^ evolution of the Wigner function and verify quantitatively that erasing of the off-diagonal elements 

corresponds to averaging out of the momentum-space details in the Wigner function. 
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I. INTRODUCTION 



Nonequilibrium Green's function (NGF) techniques [U-S] represent a powerful tool to de- 
scribe the evolution of correlated quantum many-body systems. While direct applications of 
those techniques have produced many interesting results in other fields [SHE] , within nuclear 
physics they have been mostly used for derivations, as in P [914T2]. rather than exploited 
directly [U [131 IE]- Moreover, the direct applications of NGFs have either pertained to 
uniform matter, see e.g. p^HTS] , or to the lowest-energy range of nuclear reactions [Ti] . 
By contrast, application of the static limit of Green's function theory to stationary nuclear 
states has been advanced much farther, accommodating different types of many-body cor- 
relations in the description [19]. The situation with the nonequilibrium theory in reactions 
may be, in part, attributed to the serious numerical issues faced by that theory. The im- 
mediate obstacles, however, will become less of an issue as numerical capabilities increase. 
Ultimately aiming at a direct application of the NGFs to the reactions, we consider here 
strategies for handling the challenges ahead by considering reactions of nuclear slabs in one 
dimension. The examples that we shall discuss concern the mean-field approximation of 
the NGF formalism. In the future, we shall explore the extension to correlated dynam- 
ics and higher dimensions, relevant for a realistic description of nuclear dynamics. This 
first study also provides insights into the dynamics of density matrices, coarse-graining and 
time-reversibility in quantum mechanics. 

Even with correlations incorporated into the Green's function approach, that approach 
is not going to be free of important limitations. Thus, the primary quantities within the 
Green's function theory [U [2] are single-particle functions. Given the general practical diffi- 
culties, it is not likely that an approach that relies on more-body functions as independent 
quantities can be soon developed for nuclear reactions. With this, the effects of correlations 
on the dynamics of single-particle functions may be accounted for on the average, assum- 
ing that the effects of correlations can be themselves expressed in terms of single-particle 
functions. Specifically, within the Green's function theory [HE], the single-particle functions 
satisfy Dyson-type equations in terms of single-particle self-energies. In the differential form, 
those equations are known as Kadanoff-Baym (KB) equations When the self-energies in 
those equations are approximated in terms of the single-particle functions, the equations be- 
come closed. Theoretically consistent simulations of the dynamics of correlated many-body 
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systems can, in principle, be obtained from these equations. The approximations may ac- 
count for the effects of different types of correlations, but only on an average, single-particle 
level. Not surprisingly, in the semiclassical limit the KB equations yield the Boltzmann 
equation (BE). 

With the inherent averaging over the more-body effects, the Green's function approach is 
likely to be more suitable for central nuclear reactions than for peripheral ones. In the for- 
mer kind of reactions, many particles participate and interactions between those particles are 
repeated over and over. An averaging over many-body effects takes place physically. The de- 
scription of central reactions could, thus, naturally benefit from a practical implementation 
of NGF techniques. These developments might potentially improve our understanding of 
reaction processes such as fusion and deep-inelastic collisions at low energies, multifragmen- 
tation at intermediate and high energies, and vaporization of the participant zone at high 
energies. Interestingly, the effect of fluctuations might eventually be included in the Green's 
function description by using stochastic methods pUl fIT\. 

Historically, just a handful of methods have been developed to describe central nuclear 
reactions, that simultaneously exhibited some generality and could also be employed in 
generating practical predictions. The time-dependent Hartree-Fock (TDHF) method has 
been exhaustively employed in describing low-energy reactions [22H23]. The semiclassical 
BE has been extensively used to analyze reactions at intermediate and high-energies [251 126] . 
Moreover, molecular approaches, sharing elements of both TDHF and BE, have also been 
successfully applied for various nuclear reaction purposes |2B] • 

The semiclassical BE and its related approaches have emerged as a way of dealing with the 
growth of complexity of reactions with increasing incident energy, while exploiting the short- 
ening of the particle de Broglie wavelengths. Because of their semiclassical nature, however, 
these descriptions have remained genuinely disconnected from the methods employed for 
nuclear structure, peripheral reactions or giant nuclear excitations, all with quantal under- 
pinnings [29]. More importantly, there is no systematic way of improving upon the BE- type 
approaches. Their accuracy has remained elusive, since direct comparisons to data necessar- 
ily involve adjustable parameters. At times, cross-comparisons between different approaches, 
such as molecular dynamics and TDHF, have been attempted [30] . 

TDHF [23] (and some of its extensions [201 EI]) is an approach to central collisions that 
emerges from first principles, without ad hoc assumptions put in. The nuclear system is 
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described in terms of a product wave-function, and the TDHF trajectory is found from 
a variational principle. Individual wave-functions are seen to satisfy Schrodinger-like single- 
particle equations of motion, with potential mean-fields evaluated self-consistently. Within 
TDHF, the only allowed nuclear excitations are those describable in terms of the evolution 
of single-particle wave-functions. The validity of TDHF requires a negligible role played by 
correlations in the dynamics [H]. In a fermion system, antisymmetrization of the many- 
body wave-function can suppress the correlation admixtures brought in by the interparticle 
interactions, beyond what can be absorbed into the renormalization of those interactions. 
However, as the incident energy increases, the effects of Pauli principle weaken. Correlations 
can then lead to a fast thermalization of the occupation of single-particle states and to en- 
hanced stopping, compared to what can be found in TDHF. Specifically aiming at increasing 
incident energies, it is therefore important to develop a quantal approach to central nuclear 
reactions that, besides the mean-field effects, can also incorporate correlations. 

Aside from reactions, time-dependent descriptions can lead to improvements in the un- 
derstanding of nuclear structure [23] • Thus, the response properties of nuclei can be studied 
by exciting the ground state with an external field and simulating the subsequent system 
evolution [32]. The inclusion of correlations, particularly if they give rise to dissipation, can 
introduce substantial modifications in the structure of the response. Obviously then the 
correlations can significantly affect the properties of giant resonances and of other collective 
excitations and may be, in this context, considered from a time-dependent perspective [21]. 

The study of ID nuclear systems (nuclear slabs), as undertaken here, is likely to appear 
academic. However, slabs offer an excellent initial testing ground for the ideas involved 
in the time evolution, whether for correlated or uncorrelated systems. ID calculations 
are straightforward to implement numerically with nearly no computational compromise 
and, thus, serve as an optimal starting point for the time-dependent studies. Historically, 
the TDHF approach has followed such a path and the pioneering work on nuclear slabs of 
Ref . p2] became a benchmark for the large number of later studies in 2D and 3D [33] . One- 
dimensional collisions in TDHF were subsequently discussed from the semiclassical point of 
view and the dynamics were directly compared to that in the Vlasov description [30| l3l] . 
Different attempts to introduce the effects of correlations into ID nuclear dynamics have been 
described in the literature, mostly in the context of symmetric collisions of slabs. Specific 
efforts involved extending TDHF along the lines of the relaxation-time approximation [55V 
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and further the extension motivated by NGF theory, in Ref. [38]. Moreover, in Ref. [39] 
coUisions of a slab of a finite extension with semi-infinite matter were considered. 

The NGF approach that we try to develop has a chance to generalize both the TDHF 
and BE approaches. The method can also remove some of the serious limitations of those 
two approaches. Unfortunately, even at the single-particle level, NGF techniques require 
handling vast amounts of information that can easily overwhelm the capabilities of comput- 
ing systems, rendering this approach impractical. One of our primary motivations for this 
study is finding whether all the information in the Green's function is equally important for 
the reaction dynamics. Another motivation is finding out whether separate numerical in- 
frastructure needs to be developed for preparing the initial states of nuclei and for following 
the reaction dynamics. 

In the next section, we discuss the KB equations, orienting the discussion around ID 
applications. After this, in Sec. |III[ we consider the mean-field approximation for the KB 
equations as well as the basic conceptual and numerical details of the procedure we employ 
in solving those equations. The preparation of the initial state by adiabatic switching on 



of the interactions is described in Sec. IV After the initial-state slabs have been prepared, 
the collisions of slabs may be simulated. The phenomenology of slab collisions is reviewed 
in Sec. |V} The suppression of off-diagonal elements in the density matrix is of relevance 
for the practical adaptation of NGF techniques to nuclear reactions. A suppression scheme 



and some of its consequences for slab collisions are presented in Sec. VII In addition to the 
practical consequences, we discuss the loss of time reversibility, induced by the introduced 



off-diagonal cuts, in Sec. VIII The evolution of the colliding system in terms of Wigner 



function and impact of the cuts on the Wigner function is discussed in Sec. IX Finally, 
in Sec. [X], we sum up the paper and draw some conclusions, based on the results found, for 
future developments of correlated NGF approaches to nuclear collisions. Two appendices 
are dedicated to some practical aspects of ID calculations. In Appendix |X| we show how 
to interpret ID results in the 3D manner and how to directly adapt 3D mean-field param- 
eterizations for ID calculations. In Appendix |B} we discuss the choice of frequency for the 
harmonic-oscillator potential that may be used at the start of adiabatic conversion of the 
single-particle Hamiltonian, when aiming at the mean-field ground-state. 
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II. KADANOFF-BAYM EQUATIONS 



The KB equations describe the evolution, under rather hberal assumptions, of the ex- 
pectation values of products of two single-particle annihilation and creation operators at 
different time arguments, i.e. 2-point Wightman functions. In the context of nonequilib- 
rium theory, these are related to the time-ordered single-particle Green's functions and give 
access to the most important observables of the reaction. Within a NGF picture, the initial 
state of the system is specified in terms of an A-hodj density matrix, pi, at time to. Such 
an initial state is required to be of some simplified structure. Here, we shall assume that the 
initial state is uncorrelated, i.e. completely describable in terms of one-body density ma- 
trices. This approximation can eventually be relaxed, leading to a modification of the KB 
equations, which we shall not treat here [2]. Under the mentioned constraints, the Wight- 
man functions are defined as expectation values with respect to the initial density matrix. 
Pi, of products of Heisenberg-picture creation, a^{x,t), and destruction, d{x,t), operators: 

h; X2, ^2) = i {a\x2, ^2) a{xi, ti)) , (1) 
G^{xi, h; X2, 12) = -i (a(xi, ti) a\x2, ^2)) , (2) 

where discrete quantum numbers are suppressed, fermions are assumed, and (■) = Tr {pj ■}• 
Up to a factor, the time-diagonal [i.e. t = t') Green's function reduces to the one-body 
density matrix of the system. 

The KB equations for ID systems, governing the evolution of Green's functions in their 
arguments, 

+ /" dis+(ii)G^(ii') + / dis^(ii)G-(ii'), (3) 
+ /" diG'+(ii)s^(ii') + / diG'^(ii)s-(ii'), (4) 

Jto Jto 

follow from considerations of the equations of motion for creation and destruction opera- 
tors [1]. The simplified notation 1 = (xi, ti, ai, . . .) has been introduced and the retarded 
and advanced (-) functions are defined according to: 

2) = F\l, 2) ± e [±{h - t2)] [F>(1, 2) - F<(1, 2)] , (5) 
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with standing for a possible singular contribution at ti = t2- The generalized self-energy 
S(l, 2) introduces interaction effects on the time evolution and also describes excitation pro- 
cesses within the system [H The self-energy in the previous equations has been separated 
in two different components. The first one involves the Hartree-Fock (HF) contribution, 
S//i?(l,2), which accounts for the instantaneous, one-body interaction of the considered 
particle with the mean-field produced by the other particles of the system. The term involv- 
ing describes time- dependent excitation processes, beyond the mean-field changes. Such 
terms account for the effect of correlations on the time evolution and need to be included 
for a complete description of nuclear reactions. 

While we orient our discussion around the intended ID applications, the extension of 
the Green's functions and of the KB equations to 3D is obvious. Otherwise, the complex 
integro-differential KB equations have to be solved in a self-consistent manner, since the self- 
energies are functionals of the Green's functions that are being solved for [ID]. For certain 
intrinsically-consistent many-body approximations to the self-energy, one can show that the 
time evolution induced by the self-consistent KB equations preserves the conservation laws 
obeyed by the system as a whole |2l HQ]. This is not a trivial issue, since the conservation 
laws can be broken by many-body approximations that, on the face, have quite sensible 
appearance. On the other hand, when obeyed, the conservation laws can be extremely 
useful in practice, as e.g. in testing the numerical implementation of the equations. 

An attractive feature of the KB equations is their generality. The time evolution induced 
by Eqs. Q and Q can easily incorporate various types of correlations, described by different 
approximations to the self-energy [1] . In other words, the NGF framework is able to include 
systematically more and more complicated processes in the self-energy without spoiling the 
conservation properties. The KB equations are nominally derived without any particular 
assumption on the physical system under consideration. Consequently, they have been used 
to study the time evolution of a number of many-body problems. In the nuclear context, 
only the time evolution of uniform nuclear systems has been discussed so far [U [13] . We are 
not aware of any attempt to generalize these works to finite nuclei. Elsewhere, the KB 
equations have been used to study the uniform electron gas [H] , semiconductor systems [5] , 
nanostructures [7], inhomogeneous atomic and molecular systems [8], or quantum dots |12] . 

The KB equations account for retardation effects if terms beyond the mean-field are 
included in the time-evolution of the system. From Eq. (|3]), one can easily see that the 
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Green's function at times ti and tii depends on the Green's functions and self-energies 
at all the previous times t, to < ^ < and < t < ti/ , via the time integrals on the r.h.s. 
Consequently, to find a solution of the KB equations, one must keep track of all the previous 
time-steps. This might become a major concern in the numerical implementation of the KB 
equations |13]. In the mean- field approximation, under which the effects of interactions are 
local in time, this is not a practical hindrance. 

Finally, it is important to note that the KB equations generally respect the non-local fea- 
tures of quantum mechanics. The matrix structure of the single-particle functions, in space 
and time, codes information about single-particle positions and momenta, as well as about 
energies. In nonstationary and nonuniform systems, neither energies nor momenta are well 
defined. In the self-energy terms in the equations, both the nonlocalities associated with the 
finite range of interactions are accounted for, as well as those of de Broglie type [2]. 

The matrix structure, effectively doubling the number of variables compared to the TDHF 
approach [22], puts the NGF approach at a computational disadvantage. Let us for example 
just consider the difficulties in storing density matrices rather than wavefunctions, in D 
dimensions. A uniform mesh of size Nx in each direction will yield locations for which 
wavefunction values need to be stored. With A^^ occupied single-particle states, the wave- 
function information for TDHF can be stored within a A^ x A^^ matrix. On the other hand, 
the density matrix needs to be stored in a generally much larger N^^ matrix. The difference 
in storage, between NGF and TDHF, is, though, not yet large in ID. Even for fairly accu- 
rate meshes with A^. ~ 100, the numerical implementation of NGF is not computationally 
demanding in ID and can be carried out straightforwardly. In fact, we expect that even the 
correlated calculations may be carried out in ID without any truncation compromises. 

Inclusion of correlations in NGF requires, further, augmenting the dimension of the stor- 
age matrix to account for the time variables (ti, ti') |13]. Overall, because of the high-power 
growth with the number of dimensions D, the storage becomes a serious issue for D > 1. 
The calculations at D > 1 likely need to involve truncations in space-time meshes. To opti- 
mize these truncation schemes, a better understanding of the role played and the structure 



of off-diagonal elements in the Green's functions is essential (see Sec. VI for further details). 
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III. MEAN-FIELD APPROXIMATION 



The KB equations simplify substantially when the correlation effects, described in terms 
of S^, are neglected. In that case, the evolution equations for and can be decoupled. 
Since the single-particle observables are more straightforwardly expressed in terms of G^, 
we concentrate on the evolution of the latter. The evolution equations simplify even more 
when assuming a negligible range for the nucleon-nucleon interactions. The self-energy then 
takes the form 



SH^(ll')=5(l-l')f/(l), 
where U depends on local densities. With Eq. (|6|, the KB equations for G< become 



(6) 



d 

ih—G<{xi,ti;xi>,ti, 
oil 



d 

-ih—— G<{xi,ti;xi>,ti 

Oty 



2" + U{Xii,ti>) 



G<{xi,ti;xv,tv) , (7) 
G<{xi,ti;xv,tv) . (8) 



2m dx\, 

Thanks to the instantaneous nature of T^up, the set of equations for the time-diagonal 
elements of the Wightman function, ti = tii, can also be closed. Up to a factor, those 
functions are identical to the single-particle density matrix p, 



p(xi,xi']t) = —iG'^{xi,t]Xi',t) 



(9) 



A combination of Eqs. ([T]) and ([s]) yields: 



-.9 , 

ih—p{xi,xi']t) 



r d 



+ U{xi,t) -U{xv,t) 



p{xi,xv]t), (10) 



2m 1^ dx\ dx\, 

which describes the time evolution of the density matrix in the mean-field approximation. 

The connection between the equations presented here, i.e. the mean-field approximation 
of the NGF approach, and TDHF, is established by decomposing the Wightman function in 
terms of Ng single-particle states: 



Ns 



-iG^{xi,ti;xi',ti>) = ^0a(a;i,ti)0*(xi/,ti 



(11) 



a=l 



A substitution of (11) into ([T]) then yields a set of A^^ TDHF single-particle equations: 



9 , , , 

lh—(pa{Xl,t) 

at 



2m dx\ 



+ U{xi,t) 



4>a{Xi,t) . 



(12) 
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At equal time-arguments for the Wightman function, the decomposition (11 ) reduces to that 
for the mean-field density matrix ([9]), 

p{xi, xr;t) = ^ (paixi, t) (f)*^{xr,t) . (13) 

a=l 

The equivalence between the two approaches implies that, for the same initial conditions 
and mean-field interactions, both the mean-field NGF and the TDHF equations will drive 



exactly the same time evolutions. Notably, the decomposition (11) is possible, in terms of 



a finite number of terms, with normalized to unity, when the decomposition of the density 



matrix (13) holds at some stage of the system evolution, such as for the combination of 
ground states prior to a collision, and the system otherwise follows mean-field dynamics 
that does not change the normalization of the states. 

In the following, we shall assume a spin-isospin saturated system. The Wightman func- 
tions, and the density matrix in particular, are then diagonal in the spin and isospin indices 
(which are still suppressed). Further, the diagonal values of the functions are independent of 
those indices. In the strict ID interpretation of the ID calculations, the density of nucleons 
is 

ni{x,t) = u p{x,x;t) , (14) 

where z/ = 4 represents the spin-isospin degeneracy of the system, and the nucleon number 
is A = uNg. A 3D interpretation is further possible where the matter in 3D is assumed 
uniform in two directions and nonuniform in the third x-direction. The y and z directions 
are described in terms of a set of plane-wave wavef unctions, independent of the a state 
or of time. With this, the density in the 3D interpretation, n3{x,t) = n{x,t), becomes 
proportional to the ID density, 

n{x,t) = C,ni{x,t) . (15) 

The advantage of developing a 3D interpretation is that it allows to employ well-known 3D 
mean-field parameterizations for ID calculations. 

In Appendix |X| we arrive at the following expression for the scaling factor: 



3 V 6z/2 
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by demanding that, within the 3D interpretation, the system energy minimizes at the normal 
nuclear density no- We choose a simple mean- field parametrization for our calculations: 

U=^to n{x, t) + ^h [n{x, t)]''"' , (17) 

with the parameters to = -2150.1 MeV-fm^ = 14562 MeV-fm^"*'^'" and a = 0.257 fitted to 
the saturation properties of nuclear matter: no = 0.16 fm~^, energy per nucleon of —16 MeV, 
and incompressibility of 220 MeV. 



With the potential U dependent on density in (17), the equation of motion (10) for the 
density matrix becomes nonlinear. However, the nonlinearity nominally concerns just the 
diagonal x = x' elements of the matrix. Any features of the matrix that are far away from 
the diagonal, with relatively compact support, will evolve over finite times in such a manner 
as if the equation were linear. Even for more realistic mean fields, dependent on off-diagonal 
matrix elements, the region of that mean-field dependence would be still limited to the 
immediate vicinity of the diagonal, making the equation in practice linear for features of the 
matrix far away from the diagonal. 

For a local mean-field, the time evolution of the density matrix can be implemented 
numerically in a rather straightforward way following the so-called Split Operator Method 
(SOM) |31]. Over a time-step At, the time evolution of the density matrix described by 



Eq. ( 10 ) is formally given in terms of the single-particle kinetic K and mean-field U operators: 
p(t + At) =T^{e-^r"'d*'[^^+^i]} p(t)T'^|et/^"'i*'[^^'+'^^']} , (18) 

where T'^ and T° are chronologically and antichronologically ordering operators, respectively. 
One problem in applying the evolution operators is that the kinetic-energy and mean-field 
operators do not commute. However, an expansion of the functions in the short time step 
At leads to the identity 

that helps to circumvent the issue of commutation over short time steps At. In the above 
expression, the evolution operator has been factored into two mean-field and one kinetic en- 
ergy factors. Within the spatial representation of the density matrix, the mean-field factor 
reduces to a c-factor. On the other hand, within the momentum or wavevector represen- 
tation of the matrix, the kinetic-energy factor becomes a c-factor. In each representation, 
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those factors are just phase factors and, for small At, they can be trivially applied to the 
density matrix. The only serious practical issue that remains is that of the switching be- 
tween configuration-space and wavevector representations. This switching can be optimally 
accomplished using Fast Fourier Transforms (FFT) [44j. Note that, in applying the evolu- 



tion factors from (19) in (18), the diagonal elements of the density matrix remain unaffected 
either in the configuration or wavevector space. However, the changes in the off-diagonal ele- 
ments, together with the Fourier transforms, combine to induce modifications in the diagonal 
elements as well. 

At each time step of evolution, the density matrix provides full single-particle information 
on the system, including density values, in either interpretation, ID or 3D. The total particle 
number in the ID interpretation is equal to 

A{t) = u j dxp{x,x;t) . (20) 

Upon rescaling with ^, that number provides further the number of nucleons per transverse 



area in the 3D interpretation. An evolution algorithm based on Eq. (19) ensures numerical 
conservation of the particle number. The kinetic energy may be either computed from the 
density matrix in the spatial or wavevector representation: 

m = 1; - p(-. = 1^ - /* pik, k; t) . (21) 



On the other hand, the potential energy is obtained from the density, according to 



V{t) = ^ / dx 



(22) 



The net energy, E{t) = K{t) + V{t), is conserved by the mean-field evolution equations, 
even though K and V may individually change with time. 

Our numerical code, that implements the time evolution for ID slabs following the SOM, 
employs a constant time step and an evenly spaced mesh in space. The typical step com- 
bination has been At = 0.5fm/c and Ax = 0.25 fm. For FFT, it is necessary to adopt pe- 
riodic boundary conditions at the extremes of the computational mesh. Typically, we have 
employed a mesh spanning the interval —L < x < L, where L = 25 fm. The code has 
been tested in different ways, such as in reproducing the analytic evolution of free Gaussian 
wavepackets. For any initial conditions, the code conserves particle number with an accuracy 
that is only limited by machine precision. The accuracy of energy conservation depends on 
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At and will be discussed in the context of manipulations of the evolution, see Sec. |VT| Some 
further testing of the code has involved verification of time reversibility, where the evolution 
was followed until a reaction has been largely completed and then ran back in time to reach 
practically the same initial conditions. Some information pertinent to that testing will be. 



again, presented in the context of manipulations of the evolution in Sec. |VIII[ Finally, we 
have tested the NGF code against a traditional TDHF code, arriving at consistent results, 
within the accuracy of those codes. 



IV. PREPARATION OF THE INITIAL STATE 



The initial state for a reaction needs to be constructed out of the ground states of the two 
nuclear systems entering the reaction. Observables for a ground state should, in principle, 
not change with time. However, if inconsistent approximations are employed for preparing 
the initial state of a system and for advancing its time evolution, spurious changes are likely 
to appear for observables that should have been static originally. Such spurious changes are 
likely to impair the proper theoretical understanding of the nuclear reaction dynamics. 

One way of preparing the initial system in a consistent manner within the NGF evolution 
is through an imaginary-time evolution [U [T3], employing analogous approximations for 
the self-energy within the imaginary- and real-time domains. The complication, though, 
is that the boundary conditions for the imaginary evolution tie extremal instances in the 
evolution. Because of this, the evolution has to be carried out in passes, until consistency is 
reached over a time interval that may possibly need to expand with the progress of iterations. 
As a consequence, the imaginary-time evolution requires a development of a numerical code 
separate to that for the real-time evolution, of likely greater complexity than the real- 
time code. Naturally, it may be worthwhile to explore alternative ways of preparing the 
initial state. When trying to avert the development of a separate code for finding such 
an initial state, an obvious option is that of an adiabatic switching from a precursor, simple 
Hamiltonian to the more realistic Hamiltonian of interest. Here, we explore that possibility, 
switching between an external potential and self-consistent mean-field. 

At time t — )■ — oo, whether ultimately aiming at a correlated or uncorrelated state, 
the system may be assumed to represent the ground state for an external potential Uq, 
such as a Woods-Saxon or Harmonic Oscillator (HO). Here, we choose the latter. Provided 
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that the mean-field interaction is switched on slowly enough, while the external potential is 
extinguished, the system may be evolved to the interacting ground state [IS]. For the sake 
of minimizing the changes in the system, the precursor potential Uq should be chosen such 
that the system's geometric characteristics do not evolve significantly, while the interactions 
are switched on. 

With the adiabatic switching on, the single-particle potential acquires an explicit time 
dependence, and we adopt 

Ut = F{t - To) Uo + [l- F{t - To)] U , (23) 

where U is our mean-field parametrization and, further, 

Uo = lmn^x^ (24) 

is the precursor HO potential. A choice of the frequency f2 to minimize the evolution of 
geometric characteristics is discussed in Appendix [B] The switching function F above is 
equal to 1 at some initial argument ti (< 0) and to at some final argument t/ (> 0) 
when the switching over is completed. This switching function is constructed in terms of 
a monotonically decreasing function / (such as jif) = —t in the linear case): 



for the argument values ti < t < tf. In the switching function for Eq. (23), we have most 
often employed 



within Eq. (25). Here, r represents a transition time that, for the sake of adiabaticity of the 



switching, should be longer than any characteristic time of the system. Notably, whenever 



f{ti) ~ 1 and f{tf) ~ 0, such as for \tij\ ^ r in the case of (25), then the two functions / 
and F practically coincide, 

m ^ fit) . (27) 

The slower the process of switching over is, the better an approximation to the mean- 
field ground-state the final state of the adiabatic evolution is likely to be. In the top panel 
of Fig. [!} we show the evolution of the energy per nucleon for a slab initiated at time 
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t = — lOOOfm/c with Ns = 2 HO shells filled {i.e. A = 8 nucleons in the ID interpretation). 
We consider different transition times, r, with a fixed latent time Tq = — 600fm/c. For any 
of the employed transition times r > 5fm/c, the energy evolves to a value very close to that 
for the static Hartree-Fock solution. Between the employed values of r = 5 and 40fm/c, the 
energy per nucleon obtained at t = changes just by mere 0.04 MeV! Judging the quality 
of the approximation on the basis of the energy alone, though, can be treacherous. This 
is because the energy is quadratic in the deviation of wavefunction from the ground state, 
around the energy minimum. As a consequence, final states for the adiabatic evolution 
might be found, with a wavefunction poorly approximating that of the mean-field ground- 
state, but giving an energy close to the ground-state value. On the other hand, density n is 
linear in the deviation of the wavefunction from the ground-state, so it may provide a better 
measure of the the wavefunction quality. With this in mind, in the bottom panel of Fig. [T| 
we show the evolution of the size of the slab, defined as 

D = 2{\x\) =2 j dx\x\n{x,t) I jdxn{x,t). (28) 

Indeed, for r = 5 and lOfm/c, slab sizes exhibit significant oscillations in the final state, 
indicating that the ground state has not been yet satisfactorily reached. For r > 30fm/c, 
as expected from the adiabatic theorem, the oscillations become insignificant, with the slab 
sizes practically coinciding with that for the static solution. Note also, in the context of 
Fig. [H that the emerging total thickness of the self-consistent slab, 2D ~ 6.2 fm, is quite 



close to that expected on the basis of Eq. (All), i ~ 6.4 fm. 



To supplement the above results, we show in Fig. |2] the evolution of density from the 
predecessor HO form to the self-consistent one, for the transition time r = 40fm/c. The 
densities shown here and throughout this paper will be those of the 3D interpretation. 
Already at t ~ — 400fm/c, the asymptotic form of the density is reached. Notable within 
the density is the dip at the center of the slab. This is a refiection of the node in the second 
asymmetric orbital. As such, the dip represents a shell effect that gets to be somewhat more 
pronounced for the mean-field than HO potential. The average density across the slab center 
ends up not too far from the normal density of uq = 0.16 fm~^. Overall, during the adiabatic 
switching on of the mean-field, the changes in the density are quite modest, attesting the 



utility of Eq. (B4) in choosing the HO frequency, f2. 



Let us also mention that the switching function F with a Fermi-Dirac / as shown in 
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Eq. (26) is in practice optimal. This particular switching function leads to a close approxi- 
mation of the ground state within the shortest possible time, from among tested functions. 
To illustrate that point, we show in the top panel of Fig. [3] the evolution of the energy per 
nucleon and the width for the slab considered previously, initiated with Ng = 2 HO shells 
filled at — lOOOfm/c. Different switching functions are considered: one of the Fermi-Dirac 



type as in Eq. (26), with r = 40fm/c, and, further, a linear and a piecewise quadratic 
function with suitably chosen characteristic transition times. In each case, the bulk of the 
changes occurs over the same period of ~ 400fm/c and is centered at tq = — 600fm/c. 
A first striking observation concerning Fig. [3] is how closely the energy per nucleon follows 
the shape of the switching function. This can be attributed to the relative adiabaticity of the 
switching, with the deviation of the wavefunction being small from an instantaneous ground 
state, and to the quadratic dependence of the instantaneous energy upon that deviation. 
As evident in the bottom panel of Fig. [3| for the Fermi-Dirac F also the width of the slab 
follows closely the shape of the switching function, remaining in particular stable at late 
stages of the transition. On the other hand, the width of the slab for the linear F, in the 
same panel, oscillates both during and after the transition, underscoring a poorer adiabatic- 
ity of the transition for that F. Inferior adiabaticity for the linear F may be surprising, as 
the slope of the linear / is lower than the slope of the Fermi-Dirac / at every t. Even in 
the case of the piece-wise quadratic F, oscillations in the width may be discerned for the 
final state of the transition, although they are of significantly lower amplitude than for the 
linear switching function. This points towards the importance of using smooth functions for 
adiabatic transitions. 

We hope that, in a similar manner as with the switching on of the mean- field, the retarded 
self-energies S^, representing correlations, can be switched on within the KB equations, 
with parallel changes taking place in the instantaneous self-energy. The goal would be, 
as here, for the system to stay close to the ground state for a given instant of evolution. 
With correlations, however, the switching would need to be slow compared to characteristic 
correlation times |16]. The success of the switching procedure could be, in particular, judged 
again with the degree to which the final state is stationary. 

After the density matrix for an interacting ground state has been constructed, it needs 
to be boosted in order to be incorporated within the initial state of a collision. The boost 
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is accomplished through a simple multiplication by momentum phase factors 

Pi(xi,xiO =e^''"^/Vo(a:l,a;lOe-^^^l'/^ (29) 

Here, pi is the density matrix for a slab moving at momentum per nucleon P while po is 
the matrix for an idle slab. As a consequence of the boost, the kinetic energy for the slab 
increases by the amount AK/A = P'^ /2m. The net density matrix is constructed as a sum 
of the density matrices for two countermoving slabs: 

p{Xi, Xv) = pi{Xi,Xi') + P2{Xi,Xi') . (30) 

In this paper, we present symmetric collisions, with the second slab being a symmetric 
reflection of the first. With this, the density matrix of the second slab is related to the 
matrix of the first by 

P2(a;i,a;i/) = pi(-a;i, -a;i/) . (31) 

Under the reflection symmetry, the center-of-mass (CM) energy of the collision becomes 
simply Ecm/A = P^/2m. 



The net density matrix as a sum of the two individual density matrices (30) is equivalent 



to the Slater determinant for the two nuclei being built up in an independent manner 
With this, no coherence is assumed between the initial, far-away states of the two nuclei. 



V. COLLISIONS OF SLABS 



The relatively simple ID model of nuclear collisions discussed here demonstrates a sur- 
prisingly rich range of phenomena [221 EB]- Qualitatively different physical processes are 
observed within the model when changing the CM energy for the reactions. At low collision 
energies {Ecm/A ~ 0.1 — 0.5 Me V), the nuclear slabs fuse into one compound slab that re- 
mains excited for long times. For intermediate energies {Ecm/A ~ 0.5 — 15MeV), a fusion 
process is observed, followed by a break-up into a number of smaller slabs. Higher reaction 
energies {Ecm/^ > 15 MeV) yield a pile- up of density at the system center, followed by a vi- 
olent break-up phase with the formation of a fragmented low-density neck. The process is 
reminiscent of multifragmentation in nuclear reactions. Since analyzing phenomena in slab 
collisions is not the principal goal of our paper, the discussions of those different phenomena 
will be quite limited here. 
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At low collision energies, the left and right slabs approach each other slowly until they 
get in contact and start to overlap. In our simplified model, without Coulomb effects, 
fusion can take place no matter how low the collision energy. The formed compound slab 
acquires a total mass of 2A and an excitation energy that exceeds the collision energy. 
As an example. Fig. |4] displays the evolution of the density for two A = 8 slabs colliding 
at the rather low collision energy of Ecm/A = 0.1 MeV. At t = 0, the slab centroids 
are separated by 15 fm. Following the contact at t ~ 250fm/c, a compound slab is formed. 
The energy associated with the translational motion of the slabs is converted into the energy 
of collective oscillations for the compound slab. Those oscillations are reflected in temporal 
changes of different characteristics of the system. In particular, the density proflle is observed 
to oscillate, in the later stages of the collision, around the shape of the Hartree-Fock ground 
state of the A = 16 system, shown for reference in the last panel of Fig. |4} The oscillations 
cannot easily damp out in this model [35], because their energy cannot be transferred to the 
transverse degrees of freedom and because of the constrained values of the occupations of 
the evolved single-particle states, which inhibit the transfer of energy from longer to shorter 
wavelengths. 

To provide insight into the process of slab fusion as a function of collision energy, we show 



in Fig. [Sjthe evolution of the system size, as given in Eq. (28), for different collision energies. 
Before contact, at early times, the slabs approach each other at a constant relative velocity 
equal to the slope of the size. With the slabs starting out always 15 fm apart, higher 
collision energies involve earlier contact times, signaled by a sudden change in the slope of 
the indicated system size. As the collision energy is changed, qualitatively different outcomes 
are observed for the collision flnal states. For the lowest collision energies, at times t > 
200fm/c, the system size oscillates around a central value. Boundedness of the size signals 



the formation of a compound slab of mass 2A = 16. Following Eq. (All), the expected 
size for a ground-state slab with that mass number is i/2 = 6.4 fm and this appears to be 
about the value around which the oscillations occur at late times in Fig. |5| For collision 
energies in excess of about Ecm/A = 0.45 MeV, the two slabs interpenetrate but then 
separate back into two A = 8 fragments that move apart at a lower relative velocity than 
initially, as indicated by the reduced slope of the late time size in Fig. [5j The reduced 
velocity points to internal excitations for the separating fragments. It should be noted that, 
counterintuitively, the transition between the fusion and fusion-flssion regimes is neither 
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gradual nor monotonic. Thus, while at Ecm/A = 0.2 MeV the slabs fuse in Fig. |5| the slabs 
separate at 0.3 MeV and then fuse again at 0.4 MeV. The slabs that separate at Ecm/^ = 
0.3 MeV, after one full oscillation of the system, move apart faster than either at 0.5 MeV 
or 0.6 MeV. These observations point to underlying resonant phenomena in the transfer of 
energy between translational and internal degrees of freedom of the slabs [22]. Even when 
TDHF is supplemented with a collision term in the relaxation time approximation, the late- 
time separation back into original fragments persists over a certain range of energies |37] . 

As the collision energy increases further, the final states of the collisions undergo ad- 
ditional changes. Thus, at energies Ecm/^ ^ 1.5 MeV, the system subdivides into three 
fragments. At the center of the system, a slab representing a mass of about A = 8 forms and 
stays at rest in the system CM. Besides, two residual slabs of mass of about A = 4 form, 
moving symmetrically forward and backward in the CM. An example of such a collision, 
at the energy of Ecm/^ = 4 MeV, is shown in Fig. |6} Following the fusion of the original 
slabs at t ~ 50fm/c, peaks in density, not far from wq, begin to emerge symmetrically at the 
edges of the matter at t ~ lOOfm/c. By t ~ 125fm/c, those peaks separate from the central 
fragment. At first, this remainder central region seems to fragment further. However, while 
the leading peaks in the density manage to separate from the central slab, the two more 
central peaks lack sufficient energy to separate from each other. In the end, the central 
region recontracts and oscillates, representing, in this, a highly excited state of an A ~ 8 
system. The A ~ 4 fragments remain excited too, as evidenced in the small changes with 
time of their central density. 

Generally, as the collision energy increases, the maximal density reached at the sys- 
tem center rises. The time during which the density stays substantial, on the other hand, 
decreases. The decompression of the central region becomes more and more violent for 
higher energies, thus enhancing the number of produced fragments. As an example. Fig. [7] 
shows a collision of two A = 8 slabs at Eqm/^ = 25 MeV. At the maximal compression 
point, t ~ 30fm/c, the central density exceeds by more than 70% that in normal matter. 
By t ~ 50fm/c, leading density peaks begin to form at the edges of matter, not far from 
uq, similarly to the situation at Ecm/A = 4 MeV. While these peaks continue to separate, 
the central region undergoes structural changes. In contrast to the situation at 4 MeV, 
rather than attempting to fragment into pieces at normal density, the central region under- 
goes a rapid decompression and a low-density neck region forms between the leading residual 
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fragments, see panels for t = 60 and 70fm/c in Fig. [7j The neck subsequently begins to 
recontract into separate fragments. The fate of those low-density structures is difficult to 
predict without increasing the size of the calculational box. 

VI. OFF-DIAGONAL STRUCTURE OF THE DENSITY MATRIX 

At equal arguments in a specific representation, such as the configuration or the momen- 
tum representation, the density matrix or, more generally, the function —iG^ of Eq. ([T]), 
yields the single-particle density in that representation. At different arguments, the off- 
diagonal matrix elements reflect existing correlations between magnitudes and phases of 



single-particle wavef unctions, cf. Eq. (13). As we have discussed earlier, the task of fol- 
lowing all the elements in 3D is likely to overwhelm present computer storage capabilities. 
Nevertheless, the quantities of most direct physical interest, including densities, tend to be 
associated with either diagonal values of the functions or values close to the diagonal [1]. 
Hopefully, the matrix elements which lie sufficiently far from the diagonal may be ignored and 
realistic NGF calculations may actually be implemented numerically. Note, however, that 
the temporal evolution couples different matrix elements, cf. Eqs. Q and (10). Moreover, 



off-diagonal elements in one representation are needed for the transformation to another 
representation and thus for obtaining even the diagonal elements there. Any discarding of 
the elements will have to be carefully tied to a specific representation and justified within 
that representation. Impact on other representations of interest will need to be understood. 

In the following, we shall examine the off-diagonal structure of the density matrix within 
the spatial representation for colliding ID slabs. Figure [8] illustrates what may be generally 
expected, as far as the off-diagonal structure of the density matrix is concerned. The wave- 
functions of the initial slabs are confined to within the size of the respective slab. For the 
density matrix (see the left panel of Fig. [s]), this implies a compact support limited to the 
size of the slab in all directions. In particular, a square-like region in the variables x and x' 



will develop around the x = x' axis, cf. Eq. ( 13 ) for equal time-arguments. In the late stages 



of a higher-energy collision, such as in Fig. [6| the single-particle wavefunctions are shared 
between different fragments, reflecting the fact that the individual original nucleons had 
been probabilistically distributed within the different fragments. Correspondingly, the den- 
sity matrix is likely to develop a patch-like structure extending far away from the diagonal 
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(see the right panel of Fig. |8j). As far as magnitudes are concerned, the contributions to 
the diagonal elements from different wavefunctions are all real and positive, as evidenced in 



Eq. (13). On the other hand, phase differences between wavefunction values at different lo- 
cations will make the contributions that arise from those individual wavefunctions come with 
different phases to off-diagonal elements, plausibly suppressing those elements compared to 
the diagonal. Larger differences in position and larger energies will generally increase such 
phase differences and one might expect corresponding greater suppression of off-diagonal 
elements. 

Physics beyond the mean-field dynamics (short-range correlations, particle and gamma 
decays, etc.) are hkely to introduce additional decoherence between the separating frag- 
ments, beyond that stemming only from different mean-field orbitals. Correspondingly, 
higher values for the off-diagonal elements are likely to persist more in a mean-field ap- 
proach than in any simulation with correlations. As a corollary, if one finds out that the far 
off-diagonal elements may be disregarded within the mean-field dynamics, then it should be 
even more possible to disregard such elements in more realistic approaches. 

The ID nature of the dynamics discussed here allows for a straightforward visualization 
of the density matrix in terms of intensity plots. Figures [9, 11 and 13 show such intensity 
plots for p{x,x';t) at different times t for the Eqa-i/A = 0.1, 4, and 25MeV collisions of 
A = 8 slabs. The top and bottom panels exhibit, respectively, real and imaginary parts 
of the matrix. Three sample times have been chosen to represent the initial, overlap and 
late instants within each collision. In addition. Figs. |4| [6] and [7] show values of the density 
matrix along given lines perpendicular to the x = x' line for the same sample times. All 
of the mentioned figures reflect the hermiticity and/or positive definiteness of the density 
matrix. With the matrix being hermitian, its real part is a symmetric function of its two 
spatial arguments, while its imaginary part is antisymmetric. In particular, the imaginary 
part vanishes along the x = x' diagonal. In addition, when studying symmetric collisions, 
we encounter the {x,x') — ?■ {—x, —x') symmetry for the matrix in the {x, x') plane. 

In the following, we discuss, in sequence, the structure of the density matrices in config- 
uration space for the collisions at the three energies specified previously. We start out with 
the matrix for the initial state of the low-energy slab fusion, illustrated in the left panels 
of Fig. |9} In addition, the top left panel in Fig. 10 shows a section across the t = density 
matrix, along a line perpendicular to the diagonal of the matrix. The single-particle orbitals. 
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used for the construction of the ground-state, are real and centered around ±7.5 fm for the 
two slabs. As it is apparent in Fig. |9j the compact orbitals yield compact supports for the 
individual contributions of the two slabs, to the total density matrix. These contributions 
occupy square-like regions in the {x,x') plane, centered around the points x = x' = ±7.5 fm. 
The density matrix peaks along the diagonal, x = x'. With the orbitals constructed inde- 
pendently for the two slabs, coherence patches at ggested by the right panel in 
Fig. [81 are absent in Fig. [91 At t = 0, positive peaks in the real part of the matrix along the 



diagonal, are followed by negative values farther away from the axis, see Fig. M and also 10 



For this reaction, the imaginary part is generally much smaller in magnitude than the real 
part, both at t = and later. 

If there were no boosts involved, the matrix out of real orbitals would have been, in fact, 
purely real. Moreover, if only the lowest, a = 1, purely positive orbitals were filled, the 
matrix would have been positive. The negative values of the real part of the matrix can 
therefore be attributed, in the initial state, to the filling of the second orbital, odd with 
respect to the origin of the slab. With the filling of just those two orbitals, the matrix in 
the central region of the slab can be already successfully compared to that expected in the 
ground-state of nuclear matter: 
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(32) 



where ppi is the ID Fermi momentum, cf. Eq. (13) and Appendix [Aj The function on the 
r.h.s. of Eq. (32) is, in particular, familiar from the context of slit diffraction. However, 
in contrast to slit diffraction, a band of momenta rather than positions is permitted here. 



The result of Eq. (32) predicts alternating regions of positive and negative values in the di- 
rection perpendicular to the diagonal as well as a central positive region twice as wide as 
the others. In particular, the first zeros in the off-diagonal direction are anticipated at 

Pfi 

which agrees reasonably well with what may be observed in the first panel of Fig. [TOj In the 
estimate, we have used Eqs. (A7) and (A9) with uq = 0.16 fm~^, yielding ppi/^ = 1.0 fm~^. 
Furthermore, when a slab is boosted to an average momentum per nucleon P, the ground- 



3.2 fm, 



(33) 



state density matrix gets multiplied by a phase factor, cf. Eq. (29), 
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With a real ground-state matrix po, a simple relation follows for the boosted matrix, 
Im pi{x,x') = tan [P(a; — a;')//i] Repi(a;, x'). Whenever the initial momentum per nucleon 
is low, i.e. P/h < 2/i, with i the total thickness of the slab, the imaginary part of the 
matrix will be small compared to its real part. In addition, in this limit, the real part be- 
comes practically identical to the ground-state matrix, justifying our analysis in terms of the 
ground-state matrix for a collision at Ecm = 0.1 MeV {P/h ~ 0.07 fm^^). With the result 



of Eq. (32) for nuclear matter, we can further expect in the central region of a ground-state 
slab: 
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In general, we see in the central expressions of Eqs. (35) and (36) that the off-diagonal 
oscillations in an initial density matrix, associated with its Fermi momentum pp^i, will gen- 
erally compete with the oscillations associated with the initial momentum per nucleon P. 



The approximations on the r.h.s. of Eqs. (35) and (36) are in the limit of P/h <^ 2/ 



The expectations from the r.h.s. of (35) and (36) qualitatively agree with what is displayed 



for the initial slabs in Figs. [9] and 10, over the supports associated with the slab size. Note 
that for the right-hand slab the initial momentum is negative, so its imaginary part is of 
the opposite sign, for a given location relative to the center of the support, compared to the 
matrix of the left-hand side slab. 

As time progresses, the slabs at 0.1 MeV get into contact and fuse. The density matrix of 
the compound slab acquires a compact support and the single-particle orbitals now spread 
out over the extent of the A = 16 system. The spreading of the orbitals is reflected in the fact 
that, upon fusion at t = 300 and 800fm/c in Fig. |9| the net density matrix becomes square- 
like, as for the individual slabs, but with a larger spread. Given the extremely low collision 
energy, the A = 16 slab is, obviously, only moderately excited. Accordingly, in the direction 
perpendicular to the the structure of oscillations in the density matrix becomes 

very similar to that for the initial ground state slabs, see Figs. |9] and [lOj This is consistent 
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with the fact that the oscillations are governed by the Fermi momentum ppi of normal 
matter and qualitatively described by Eq. (32). In fact, for the A = 16 slab, additional 
oscillations, also controlled by ppi, are seen in the direction perpendicular to x = x', beyond 
those observed for the initial A = 8 slabs. The persistence of the imaginary part of the 
matrix, at late stages, may be understood in terms of collective motion with the matter 
flowing back and forth. In fact, a finite net local flux j{x) is conditioned on the imaginary 
part of the density matrix being finite next to the x = x' axis: 
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(37) 



When moving along the x = x' axis, up to four regions are observed with alternating sign 
of the flux, see the bottom t = 300fm/c panel in Fig. |9} 

As the collision energy increases to Eqm/A = 4MeV, some qualitative changes are ob- 



served both for the initial and later-stage density matrices in Figs. [TT] and [T2| compared 
to the lowest collision energies. To begin with, the initial momentum, at P/h = 0.44 fm~^, 
begins now to compete with the Fermi momentum, ppi/h = 1.00 fm~^, in controlling the 
oscillations of the initial density matrix as a function of {x — x'), cf. Eqs. (34) and (35). 



In fact, the modulation associated with P basically eliminates the regions with negative 



values for the real part of the initial matrix in Figs. [TT] and 12 In addition, the significant 
P- value enhances the imaginary values for the density matrix compared to the previous case. 



cf. Eqs. (I34j) and (136|). 

At this energy, once the compound slab is formed, it remains in contact for a too short 
time for the single-particle orbitals to spread out evenly over the size of the compound 
system. In consequence, the square-like support for the density matrix of the compound 
system, such as at low energy, fails to develop. On the other hand, in contrast to the low 
energies, the system subsequently disintegrates into three fragments. Each of these fragment 
separately achieves a final stationary state, with approximately square-hke support around 



the diagonal, see right panels of Fig. [TT} In addition, cross-correlations develop between 
individual fragments, as signaled by the patches of significant values of the density matrix 
far away from the x = x' diagonal. As an example, the structures at a; ~ and x' ~ 15 fm for 



t = 200 fm/c in Fig. 11 , both for the real and imaginary parts of the matrix, signal the overlap 
of single-particle states between the central and the outgoing fragments. Fainter correlation 
patches may be observed between the two fast fragments, see the region at x ~ — 15fm and 
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x' ~ 15 fm there. A further insight into the features of the density matrix at the late stages 



of the colhsion, away from the x = x' axis, is provided in Fig. 12, where regions of significant 
values far away from the axis are observed, particularly at times t = (150-200) fm/c. 

The discussed correlation patches may be understood in terms of the fragmentation of 
single-particle orbitals, which we have already mentioned in the context of Fig. |8| During 
the breakup, the nucleons from the original orbitals have finite probabilities of ending up 
in different fragments. The amplitudes for those possibilities maintain a phase relationship 
leading to correlations in the density matrix. The entanglement of the internal wavefunctions 
is the sole reason for the persistence of the far-away off-diagonal patches in the density 
matrix. The real and imaginary parts for those structures are comparable in magnitude, see 



Figs. H and 12, This points to the involvement of significant relative phases, as expected 
from the significant difference in momentum and position. Note that these entanglement 



correlations persist for fragments that are 30 fm apart for the t = 200 fm/c panels of Figs. 11 
and 



In the collision at Ecm/^ = 25MeV, the initial momentum, P/h = l.lOfm ^, is already 
similar to the Fermi momentum. With this, the collision can serve as an illustration of how 



the initial density matrix is controlled by two similar momenta, see Eq. (34) and central 



expressions of Eqs. (35) and (36). The initial momentum generally controls the relative 
magnitude of the real and imaginary parts of the matrix. When the initial momentum is 



large, the imaginary and real parts become comparable, cf. Figs. 13 and 14 In addition, 
when the initial momentum becomes comparable or exceeds the Fermi momentum, it begins 



to control the oscillations of the matrix in (x — x'). In fact, in Figs. 13 and 14, the matrix 
values oscillate about twice as often as at low energies, in the direction perpendicular to 
the X = x' axis, because of the initial and Fermi momentum being about the same. On the 
other hand, from the two momenta, it is the Fermi momentum that controls the extent 
of the peak of large magnitude of the matrix around the x = x' axis, compare the left 



panels of Fig. 13 and |9j As time progresses within the Ecm/A = 25MeV colhsion, after 
the slabs come into contact, they again remain in this contact for a time that is too short 
for the square-like support to develop in the compound system. With further progress of 
the reaction, qualitatively new developments take place at this energy, not only for the 
density, but also for the off-diagonal structure. Thus, as the leading fragments recede and 
a low-density region emerges in-between those fragments, correlations develop all across the 
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calculational region, involving regions over 30 fm x 30 fm in size at t = 80fm/c. Clearly, the 
single-particle orbitals spread out quite substantially across the calculational region, with 
few spikes developing, such as those representing the leading fragments. The spikes in the 
wavefunctions give rise to different structures within the off-diagonal elements of the density 
matrix, such as ridges. 

Otherwise, for the multiply fragmented system, the general expectation is that Hubble- 
like correlations develop between the position and the average momentum per nucleon 
around that position. Consistently with such an expectation, for an orbital expanding 
in the vacuum [35], the oscillations in the direction perpendicular to the diagonal of the 
density matrix are governed, at late times, by the factor 

pix,x') oc exp y—j^) > (38) 



where and X = {x+x')/2. Comparing the above factor to that in Eq. ([34|), we see 



that momentum P is replaced by mX/t. The expectation from (38) is that of a pattern of 
hyperbolas for the zeros of the real and imaginary parts of the density matrix, corresponding 

2 /2 

to constant values of the product Xxr = ^ ~2 ■ More generally, the oscillations of the 
matrix will also follow a hyperbolic pattern. Such a pattern indeed begins to emerge in 



the t = 80fm/c panels of Fig. 13 superimposed onto the structures associated with specific 
fragments. The asymptotes for those hyperbolas are the two diagonals, x = x' and x = —x'. 



The late-time profiles of the matrix, exhibited in the last two panels in Fig. 14, are also 



consistent with the expectations based on Eq. (38). 

To summarize, the significant values of the far off-diagonal elements of the density matrix 
are associated with phase coherence within far spread-out and fragmented orbitals. Such 
a phase relationship might be of importance if the fragmented slabs were ever to recombine. 
However, at the late stages of a nuclear collision, the system tends to undergo expansion 
and progresses through a series of decays rather than through recombinations. On the 
other hand, the phase relations may be important within individual emerging fragments, 
with nucleons moving back and forth between the fragment boundaries. However, even in 
this case the correlations within the density matrix are characterized by a specific range. 



of ~ 3.2 fm, cf. Eq. (33). Monitoring of the far off-diagonal elements of a density matrix, 
substantially beyond this natural range, might be unacceptably costly within a true 3D 
system. In the following section, we examine the degree to which far-away elements can be 
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disregarded, by suppressing such elements in ID and observing the impact on the observables 
for the slab reactions. 



VII. SUPPRESSION OF FAR OFF-DIAGONAL STRUCTURE 

As we have discussed, the significant values of far off-diagonal elements of a density matrix 
can be interpreted in terms of highly delocalized single-particle orbitals. In some physical 
situations, e.g. those that are relevant for nuclear reactions, such off-diagonal elements might 
not be important. Consequently, one could expect to be able to disregard such elements 
without affecting core features of the reaction dynamics. We are going to test this here 
in practice. There are two major potential benefits of suppressing off-diagonal elements. 
The first is that the suppression procedure can provide an understanding of the practical 
role of the far-away elements in the dynamics. The second is that the procedure can introduce 
significant savings at the numerical level for the anticipated effort in realistic calculations of 
reactions. With the inclusion of correlations, and without any savings, a NGF calculation 
in D dimensions would actually require following a self-consistent evolution of N^^ x A^^^ 
values of the Green's functions. Here, A^^^ and Nt represent the number of discretization 
steps in space and time, respectively. The ability of reducing the number of matrix elements 
of the relative position, down to a nominal fraction within each dimension, is critical to 
carry through realistic NGF calculations. 

In the correlated case, it would be additionally beneficial if it were possible to consider 
only short relative times, tr = t — t', in the functions, i.e. minimize the memory timescales. 
In fact, with inclusion of relative time in the consideration, the situation gets generally 
more complicated for the spatial arguments, than in the mean-field case with tr = only. 



Thus, for tr ^ 0, a superposition of the products of orbitals, such as in Eq. (11), tends not 
to maximize exclusively around the relative position ) = 0, but rather over 

a trajectory within {xr,tr) space, that ends up at a;,. = when tr = 0. Depending on how 
strong the effect of correlations is, the values for the superposition may persist as a function 
of tr along that trajectory. Consistently, the cases of tr = and tr ^ can be handled 
within a Wigner representation for the Wightman functions, which will be discussed in this 



paper in Section IX On the other hand, consideration of short memory times only [13 



may be warranted by physically short transition times through binary interaction regions. 
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eliminating, in practice, the need to cope with the mentioned trajectories in relative space- 
time. 

In a nuclear reaction, fragments separate from each other in configuration space. It is 
therefore natural to seek a suppression of the far-away elements in the configuration rather 
than any other space where the density matrix could be represented. (Note that the density 
matrix would maximize around the diagonal in other representations as well, e.g. around 
the p = p' momentum axis). As physics observables tend to be obtained from elements of 
the density matrix close to the diagonal, any satisfactory procedure of suppressing far-away 
elements needs to retain the near-diagonal region nearly intact. Conservation laws, for one, 
rely on the near-diagonal region and can be used for testing of the suppression procedure. 
Otherwise, the results near or on the axis may be directly compared to the dynamics with 
and without the suppression of the far-away elements. 

Since only elements away from the axis need to be suppressed, x x', the suppression 
procedure needs to be in some way nonlocal. The object on which the operation of sup- 
pression needs to be carried out, i.e. the density matrix, is an operator itself rather than 
a wavefunction. Operators that act on operators are termed superoperators within the 
quantum-mechanical context [50]. These have proven useful in the studies of decoherence 
and quantum transport phenomena [51]. Thus, for example we can represent the equation 



of motion (10) for the density matrix as 



d 

ih — p = Cp, (39) 

where £ is the Liouvillian superoperator, 

C = ki + Ui- kv - Uv . (40) 

To achieve our goal of suppressing the off-diagonal elements, we supplement the Liouvillian 
superoperator with an absorptive superpotential W: 

C > Cs = ki + JJi + W{xi, Xv]t) - ky - Uy - W*{xv,Xi] t) . (41) 

The dependence on both spatial arguments in W is needed in order to suppress exclusively 
those elements that are away from the diagonal. The particular combination of the super- 
potential and its complex conjugate is enforced by the requirement that the density matrix 
stay hermitian. With the modified Liouvillian, the evolution of the density matrix over 
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a time step may be represented as 

p(t + At) =T^|e-^/^"'i*'[^^+^^+^]} p(t)T'^|ei^^"'i*'[^i'+^i'+^1} . (42) 



Elementary considerations strictly restrict the possible forms of W. From Eq. (41), 
meaningful W^s need to have real and imaginary parts that are, respectively, odd and even 
under the exchange of spatial arguments. It is also natural to impose on W the symmetries 
of the underlying system dynamics. Thus, if the dynamics is translationally invariant, 
then W should depend on the difference of arguments only, {xi — xr). With invariance 
under inversion, the superpotential needs to be an even function of (xi — xi') as well. 
Consequently, the real part of the superpotential has to vanish. To suppress off-diagonal 
elements, the imaginary part of W needs to be nonpositive. Finally, the conservation of 
particle number (or probability in other contexts) is enforced if the imaginary part vanishes 
along the Xi = xy diagonal. The simplest form of Im W , which would be a quadratic one, 
leads to a Gaussian suppression of the off-diagonal elements. That form has frequently 
appeared in the theory of decoherence [521 [53] . Other forms of IV-type functions have been 
arrived at in studying a subsystem coupled to a chaotic environment 

With our goals in mind, we chose not to suppress elements close to the diagonal for 
the matrix and to suppress strongly those sufficiently far away. We have further refrained 
from changing the suppression abruptly with location of the matrix element, to avoid any 
significant diffractive effects induced in the matrix. We have adopted piecewise parabolic 
changes for the potential, in-between the diagonal and the suppressed regions, through the 
parameterization 

, < |Xr-| < Xo , 

-iWQ2[\xr\ - XqY / dl , xo < < a;o + rfo/2, ^ ^ 

AWqIa - 2[\xr\ - {xq + do)] /dly xo + rfo/2 < |xr| < xo + rfo, 
—iW^ , Xq + dQ < \xr\ < L . 

In employing FFT within the SOM for our system, we had to make our system periodic. To 
retain the periodicity for the evolved matrix elements, we further needed to enforce period- 
icity on the superpotential. Thus, with the periodicity of our system being 2L, the super- 
potential has also been built, for \xr\ > L, to satisfy W{xr + 2L) = W{xr). 

The superpotential and the associated suppression pattern of the off-diagonal elements 



are illustrated in Fig. [15] Note that the periodicity enforces the lack of suppression for the 
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elements next to the corners at {—L,L) and {L,—L). If elements were suppressed there, 
fragments moving towards the edges of the calculational region would encounter stronger 
suppression of off-diagonal elements than fragments in the middle of the region. Effec- 
tively, the fragments would be interacting with the edges, leading to a break of translational 
invariance. 

The three parameters Wq, xq and do control the strength and shape of the superpotential. 
Specifically, Xq regulates the size of the band around the x — x' axis where matrix elements 
stay intact. The fraction of elements that get suppressed within our {x,x') plane is x = 
1 — Xq/L. The steepness in the change of W with position is regulated by do- A too small 
vahic of do niay induce spurious numerical oscillations in the density matrix, beyond those 
already inherent in the matrix. Finally, Wq regulates the pace at which the elements get 
suppressed. Over a time interval of At, the suppression factor for elements away from the 
diagonal, at \xr\ > Xo-l-do, is 6"^^°^*/^. Unless otherwise specified, we employ a rather large, 
on nuclear scale, strength l^o — 1000 MeV and a moderate — 2 fm. We have confirmed 
that, within a considerable range of values of Wq and d^, the results of the calculations are 
fairly similar. 

In the reminder of this Section, we shall attempt to quantify the practical importance of 
the far-away off-diagonal elements of the density matrix within the mean-field evolution of 
ID slabs. To this end, we shall reexamine the three coUisions discussed before, representing 
fusion, break-up and fragmentation. We shall follow the system evolution, with a progres- 
sively narrower range of retained matrix elements, i.e. smaller xq, while focussing on the 
local density and other characteristics of the system associated with the region close to the 
X = x' diagonal. Our basic interest is in how far we may be able to push, if at all, the element 
suppression without significantly affecting the region close to the diagonal. 

Before we proceed, it may be useful to cast the differential equation we solve, 

ih — p=(t + 2ilmW^p, (44) 

in one more integral form, for the sake of understanding the nature of the forthcoming 
solutions: 

dt' ei (*-*') (tp) (a;, x'- f) + e'i ^^^i^^) p{x, x'- h) . (45) 

In the above, we have made use of our conclusions on W to separate formally the time 
evolution driven by C from that associated to ImVF. The time ti precedes t. When we do 
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not suppress any elements, or else in the region where W happens to vanish, the exponential 



factors containing W become equal to 1. Otherwise, Eq. (45) demonstrates that large values 
of W act to quickly eliminate the off-diagonal elements of p. On the other hand, if p was 
not going to have significant values in the first place, where W was acting, the presence 
of W would not change the situation much for the density matrix. Finally, for any p and at 
least for finite times, in the regions that are away from the suppression, the density matrix 
is going to evolve in a manner consistent with the absence of any suppression, i. e. following 
the evolution driven by C. 

Let us start the practical investigations of the suppression with the low-energy fusion- 
reaction at Ecm/A = 0.1 MeV and consider global quantities for that reaction. The top 



panel of Fig. 16 shows the net energy of the system as well as the different contributions 
to this energy, as a function of time, for different values of the cutoff xq- The value of 
Xo = = 25fm represents the standard evolution, without any suppression. As is apparent 
in Fig. [16} the net energy remains very well conserved throughout the evolution, even for the 
relatively low cutoff value of Xq = 10 fm. Moreover, the contributions to the energy do not 
appear to depend on the cutoff for xq > 15 fm. As to the xq = 10 fm case, the kinetic and 
potential energies evolve in the same manner as for the original evolution up to t ~ 370 fm/ c. 



Thereafter, though, differences start to emerge. The bottom panel of Fig. 16 shows the 
evolution of the extent of the system for different Xq. As with the energy breakdown, the 
size appears independent of the cutoff, when xq > 15 fm. In the Xq = 10 fm case, like 
it has been discussed for the energy, the size evolves in the same manner as without any 
element suppression, until t ~ 370 fm/ c and only then differences emerge. In the lower panel, 
we also note that the differences begin to emerge only as the system begins to recontract 
following its initial expansion. In essence, the excessive, for this energy, suppression of matrix 
elements reduces the period of collective oscillations for the system. As time progresses, 
the accumulated difference in phase of the collective oscillations between the original and 
the Xq = 10 fm evolutions leads to a beating for the difference in size between the two 
evolutions. 



Figure 17 next shows the density in the Ecm/A = 0.1 MeV collision, at sample times, 
for different values of the cutoff parameter xq. In this more differential representation of 
the collision, again no differences can be observed between the standard evolution and the 
evolutions with xq > 15 fm. For xq = 10 fm, the central dip in the density is slightly filled 
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at t = 300fm/c, compared to the standard evolution. The difference at 500fm/c, between 
the Xq = 10 fm and the standard evolution, may be attributed to the different phase in 
collective oscillation. At t = 800fm/c, the densities for those two evolutions are similar as 
the oscillations within those evolutions progress through the analogous intermediate state, 



cf. Fig. 16 



Overall, we find that we can suppress up to about 50% of all matrix elements for our 
calculational region, cf. also Fig. [9} without affecting any essential features of the fusion 
reaction. When we suppress more, the system appears to be affected quantitatively at the 
later stages of its evolution. However, the system still evolves through similar stages as 
without element suppression, it just arrives at those stages at modified times. Interestingly, 
the conservation of net energy survives to a surprising degree, even when quite aggressive 
suppressions of the off-diagonal elements are introduced into the evolution. 

To an extent, we may find the above results quite understandable. If we look at Fig. |9j 
we observe that significant matrix elements extend from the diagonal only by |a; — a;'| about 



equal to the thickness of an A = 16 system, or ~ 13 fm after Eq. (All) (cf. also Fig. 16 



accounting that the thickness of a uniform system is equal to 4(|x|)). Thus, if we choose 
a large cutoff xq > 15 fm, we never suppress any significant values of matrix elements. 
A more abrupt suppression might give rise to stronger kinetic energy effects, even for small 
suppressed elements, but we have prevented that through the graduation of the suppression, 
i. e. through a finite (Iq. In the case of Xq = 10 fm, we do not suppress any significant elements 
before the slab make contact either, cf. Fig. |9j It is only after the slabs fuse and the orbitals 
manage to spread themselves over the extent of the compound slab, by t ~ 300fm/c, that 
the Xq = 10 fm suppression really begins to affect the matrix. However, phase correlations 
should only start to matter upon the return of the waves incident on the opposing boundary 
of the compound slab, hence the delay, until the system recontracts, for significant changes 
in the global features to emerge. 



Figure [18] provides additional insight into the off-diagonal structure of the matrix in 
the context of element suppression. The figure displays values of the matrix along lines 
perpendicular to the x = x' diagonal, at different xq. For the display, we have chosen the 
initial time t = and the cut at X = {x + x')/2 = 7.5 fm, across the center of the right 
slab, and, further, two later times and a cut across the system center at X = 0. With 
the t = density matrix extending up to about \xr\ ~ 7fm, neither of the suppressions in 
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the figure affects that matrix. At t = 300fm/c, the matrix has spread out in transverse 
directions up to \xr\ ~ 10 fm. The Xq = 10 fm cut starts barely to affect the matrix, while 
the Xq > 15 fm cuts leave the matrix intact. At t = 500fm/c, the undisturbed matrix is 
spread out a bit farther than at t = 300fm/c. The Xq > 15 fm cuts still leave the matrix 
intact. On the other hand, the effect of matrix suppression for xq = 10 fm has by now 
propagated from \xr\ ~ 10 fm down to the vicinity of the axis Xr ~ 0, affecting in particular 
the density. Remarkably, though, the overall structure of the matrix at \xr\ ^ Xq remains 
pretty unchanged by that last significant cut. 

In the preceding Section we have seen that the support of the density matrix spreads 
out over the whole calculational (x, x') area at higher energies, rather than staying com- 
pact as in the low-energy fusion-reactions. Given our low-energy experience in this Section, 
we might then expect that even modest suppressions of the elements in the higher-energy 
reactions could produce noticeable changes late in the system evolution compared to the 



standard evolution. To test this. Fig. 19 shows the evolution of the energy components and 
of the system extent for the reaction at Ecm/A = 4MeV, when utilizing different cutoff 
parameters Xq in the suppression of matrix elements. Consistently with our earlier experi- 
ence, we may note an excellent energy conservation in the presence of element suppression. 
However, we may also note that the energy breakdown and the system extent turn out to be 
fairly independent of the cuts applied to the elements, contradicting the naive expectations. 
Only in the case of Xq = 10 fm, we see small < 5% changes in the energy breakdown and 
system extent in the late t > 150fm/c stages of the reaction, compared to the unperturbed 
evolution. We expected that, with the increase in collision energies, the global features of 



the reactions would be more sensitive to cuts. Instead, Fig. [19] shows that these features are 
far less sensitive to the off-diagonal suppression. 



Figure 20 shows the system density in the Ecm/^ = 4MeV reaction, at sample times, for 



different cutoff parameters Xq. As has been discussed in Sec. VI, at this particular collision 
energy the system breaks into three fragments, one central with mass A ~ 8, and two moving 



forward and backward in the CM. We can observe in Fig. 20 that the cuts with xq > 15 fm 
have no practical effect on the evolution of the density. The xq = 10 fm does not modify the 
evolution throughout the formation and breakup of the compound slab. Only at later stages 
of the reaction {t > 150fm/c), slight differences emerge compared to the standard evolution. 
Thus, at t = 200fm/c, the leading fragments are ahead by ~ 1 fm compared to the standard 
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evolution. This indicates that those fragments emerge at shghtly higher velocity from the 
central region, for Xq = 10 fm. On the other hand, the central fragment at t = 200fm/c, for 
Xq = 10 fm, appears shifted in its phase of oscillation compared to the standard evolution. 
This seems similar to the situation with the fused slab at Ecm/A = 0.1 MeV. One might, 
however, wonder how an xq = 10 fm cut can affect an A ~ 8 slab, of much smaller size. After 
all, such a cut has no noticeable effect on the ground-state initial A = 8 slabs. The essential 
issue here is that the central A ^ 8 fragment is actually reconstituted from a highly expanded 



structure, cf. central panel in Fig. 20, that manages to expand perpendicular to the x = x' 
axis before recontraction, and it is therefore sensitive to the suppression in the off-diagonal 
direction. 

Additional insight into the element suppression and its consequences at Ecm/A = 4 MeV 



is provided by Fig. 20, analogous to Fig. 17 for this collision. At the earliest of the times rep- 



resented in Fig. 20 , t = 75fm/c, the single-particle orbitals and, correspondingly, the density 
matrix in the transverse direction are relatively compact. Only the Xq = 10 fm cut affects 
the matrix in a significant manner. However, the effect of the suppression, if any, would take 
a while to propagate to the vicinity of the x = x' axis. The panel for t = 125fm/c shows 
a cut along the line X = {x + x')/2 = 5.46 fm, traversing the patches within the density 
matrix that represent a correlation between the central fragment and the leading fragment 
moving to the right. The superpotentials with xq > 15 fm have had no practical effect on 
those patches, but the potential with xq = 10 fm begins to absorb them. The panel for 
t = 200fm/c, representing the cut at X = 8.74 fm, illustrates the ensuing development for 
those patches. The xq = 20 fm superpotential has had, by that time, still no impact on the 
patches, while the xo = 15 fm superpotential begins to absorb them. Finally, the xo = 10 fm 
superpotential has, at this point, completely absorbed the off-diagonal patches. The sup- 
pression of the patches, however, appears to be of no consequence for the development of 



the essential aspects of the system, cf. Fig. 19 This can be understood in the following 
terms. When the system expands, the orbitals fragment into individual pieces which take off 
with different velocities, forming fragments that move away from each other. With the frag- 
ments separating, the correlation patches in the density matrix move away from the x = x' 
diagonal. If those patches never move back to the vicinity of the axis, where they could 
contribute to the observables (momentum will be addressed later in the paper), it becomes 
irrelevant whether they are allowed to propagate forever or get suppressed. If the region 



34 



far-away from the diagonal becomes the region of 'no-return', it ceases to matter whether 
the dynamics in this region is followed or not, in order to find the essential observables. 

A similar situation is found at even higher energies when we examine the impact of el- 
ement suppression in the fragmentation reaction at Ecm/A = 25MeV. For this particular 
reaction, besides the results analogous to those for lower-energy reactions, we shall also show 
the consequences of the suppression with xq = 5 fm, that directly affects 80% of the matrix 
elements within our calculational area. This specific case will give, in particular, an op- 



portunity to discuss some practical issues involved in the element suppression. Figure 22 



shows first the evolution of the energy components and of the system extent for the different 
suppression cuts. No dependence on the parameter xq can be seen for either quantity if 
Xo > 10 fm. For xq = 5fm, only a slight deviation from the standard evolution may be seen 
within the energy breakdown at t > 80fm/c and virtually no deviation is found within the 
system extent. 



Figure 23 , with sample densities, provides more insight into the role of element suppres- 
sion at Eqm/^ = 25MeV. It is apparent in the figure that superpotentials with xq > 10 fm 
have little practical effect on the evolution of density. The likely reason is that, at this 
energy and within the followed evolution, the system rapidly breaks up into small fragments 
and no stable structures form, that might be extended in the (x — x') direction. On the other 



hand, it is evident in Fig. 23, that the superpotential with xq = 5fm affects the evolution 



of density in a nontrivial manner. At the late stage (t = 80fm/c in Fig. 23), the Xq = 5fm 
system appears to fragment into 3 rather than 5 fragments as in the evolution without any 
suppression. Some differences in the evolutions, precursing this outcome, can be already seen 



in the central panel of Fig. 23, for t = 50fm/c. Around maximal compression at t = 30 fm, 
the central dip, present in the standard density, is filled for the density with the xq = 5 fm 
suppression. A new aspect of the changes in the density with the strong element suppression 
may also be seen in the t = 30fm/c, and, to some extent, in the t = 50fm/c panel. For 
Xo = 5 fm, small-amplitude oscillations extend away from the central fused slab and take the 
density to negative values. These illustrate the fact the density matrix, a positive-definite 
operator, generally loses this property after element suppression. Indeed, while unitarity is 
preserved, i.e. the integral along the diagonal does not change, the individual values of the 
density are modified and there is no restriction on their values. We may encounter negative 
density values for less drastic suppressions as well, but these are usually quantitatively neg- 
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ligible. Generally, the appearance of those values will depend on the shape of the element 
suppression. Note that for drastic suppressions, such as Xq = 5fm discussed here, the low- 
density oscillations may become visible, even though they still make negligible contributions 
to the global quantities such as the system extent or breakdown of energy into different com- 
ponents. When evaluating potential energy, we have, so far, adopted a convention that we 
replace the argument of the energy by zero when the density is negative. In relation to the 
above, we have also found that the population of single-particle states differs from the mean- 
field values once the off-diagonal suppression is introduced. In place of the usual A^^^ = 1 
and values, fractional occupations emerge, with a distribution that could be interpreted in 
terms of a temperature. The emergence of fractional occupations can be expected, since the 
element suppression effectively decouples the portions of single-particle wavefunctions that 
have moved out with different fragments. However, for very strong suppressions, such as 
Xo = 5 fm here, we can also encounter small negative occupations as well as occupations 
slightly exceeding 1. 



Figure 25 provides next an insight into the impact of the cuts on the real part of the 
density matrix within the {x,x') plane in the Ecm/A = 25MeV reaction at t = 80fm/c. 
The matrix without any element suppression is similarly represented in the right panel of 



Fig. 13 While the superpotentials with xq > 10 fm remove matrix elements quite effectively 
there where they are applied away from the diagonal, they leave the vicinity of the x = x' 
axis practically unchanged. Even the Xq = 5 fm superpotential leaves the vicinity of the axis 



fairly similar to that of the standard evolution. Furthermore, Fig. 25 shows values of the real 
part of the density matrix, along the lines perpendicular to the x = x' diagonal, at sample 
times in the Ecm/A = 25MeV reaction, for different suppression parameters xq. As for 
the previously investigated reactions, the cuts play a more prominent role later, rather than 
earlier in the reaction. While the Xq > 10 fm cuts suppress the density matrix at the later 
stages of the reaction, they leave, as seen here again, the immediate vicinity of the x = x' 
axis unchanged. On the other hand, the xq = 5fm cut eventually induces changes in the 
structure of the matrix in the immediate vicinity of the axis, specifically at t = 80fm/c. 
Note that the different frequency of oscillations may be interpreted in terms of a different 
momentum per nucleon for a given piece of matter at a specific location. 

Overall, we have seen that extensive suppressions of far-away elements of the density ma- 
trix may be carried out without affecting central features of the reaction dynamics. The de- 
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gree to which the suppressions may be apphed depends on the energy of the reaction and 
on the size of stable structures that are formed in the dynamical evolution. The higher the 
energy and the smaller the objects, the more compact the region around the diagonal of 
the density matrix that needs to be retained. While entanglement of internal wavefunctions 
between fragments emerging from a reaction leads to correlation patches that move away 
from the diagonal (potentially up to infinity for infinite times), one does not need to follow 
the evolution of these patches if the fragments are never to meet again. As an additional test 
for the importance of the correlations, besides what has been discussed, we have followed the 
development of reactions allowing the fragments to leave the periodic computational interval. 
The fragment that leaves the computational interval enters the next one and/or enters the 
original interval back from the other side, depending on the interpretation. With this, we 
allowed the correlated fragments to collide. Within the {x, x') plane, the correlation patch 
may be seen as traveling to the diagonal within the next computational box and/or as reen- 
tering the original box from another side. On the diagonal, the patches encounter the parent 
density pulses undergoing a collision. We have carried those calculations with and without 
element suppression. In such calculations, any level of element suppression eliminates the 
patches that travel from one diagonal to another. Contrary to naive expectations, though, 
we did not find any significant differences in the evolution of those overextended collisions, 
whether or not we had suppressed the elements and the correlation patches traveling across 
the diagonals. Quite likely, the relative insensitivity to the interference in the collisions of 
the emerging fragments can be associated with the large relative momenta between frag- 
ments. There is, though, one aspect of the calculations that we found to be highly sensitive 
to any suppression of the structures in far-away elements. This is the time-reversibility of 
the dynamics, and we shall discuss that aspect of the dynamics in the following section. 



VIII. TIME-REVERSIBILITY OF THE DYNAMICS 

Generally, the time evolution of the density matrix in the mean-field approximation obeys 
time- reversal symmetry in the microscopic sense [21] , i- e. the evolution of a system, follow- 
ing the mean-field dynamics, is describable in terms of the mean-field dynamics when ran 
backwards in time. Within the symmetry, if we evolved a system past a collision, following 
the mean-field dynamics, we should be able to reverse the direction of evolution and reach 
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the initial state. We have, hke others ^9], exploited that symmetry to test our numerical 
code, and found the symmetry to be very well satisfied at the numerical level. Erasing of 
elements in the density matrix is, though, a generally irreversible process. However, the 
degree to which this irreversibility is important clearly depends on the circumstances. Thus, 
if the erased elements are small (or even zero) to start with, the erasing should presumably 
not be very important. Further, we have seen that even removal of significant elements may 
not affect much a system under expansion. In the following, we examine in practice the 
level of irreversibility induced by suppressing the elements of the density matrix, taking as 
an example the case of the parameter xq = 10 fm. This particular suppression has left the 
central features of the Ecm = 25 MeV reaction unchanged and has affected the lower energy 
reactions at a semi-quantitative but not qualitative level. 

Handling of the imaginary superpotential in the time-reversed case requires some elemen- 
tary care. At the formal level, the change of the direction of time in the equation of motion. 



Eq. (44), may be accomplished by taking a complex conjugate of the sides of that equation. 
For the conjugate equation, the forward evolution in time will be progressing through stages 
in reverse order to those in the original forward evolution. In the reversed evolution, the 
states will be described in terms of the complex conjugate matrix. In taking the complex 
conjugate of the equation, the superpotential will be conjugated as well and, as a purely 
imaginary object, it will change its sign. With this change in sign, the superpotential will 
lead to an artificial increase in the magnitude of matrix elements found beyond the cutoff, 
rather than to a decrease. However, since we are aiming at the reduction of the calculational 
space, we are interested in erasing the information no matter in which direction the evolu- 
tion progresses. Consequently, we need to change the sign of the superpotential by hand 
when reversing the time evolution. Note, besides, that a strong potential which increases 
off-diagonal elements arbitrarily would give rise to a numerically unstable evolution. 



As a characterization of the system development under time reversal, we show in Fig. 26 



the evolution of the system extent for the three collision energies considered before. Without 
element suppression, the extent for a reversed evolution retraces the extent for the system 



evolving in the forward direction. In Fig. [26| the forward and backward evolution without 
element suppression are accounted for by one line. For the forward evolution with element 



suppression, at the two lower energies, we can recognize in Fig. 26 the late-time deviations 
from the standard evolution that were already discussed in the previous section. At the 
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highest of the represented energies, no late-time deviations are actually visible. For each 
of the collision energies, the reversal of the evolution is applied at the end of the displayed 
evolutions: at t = 1000, 200 and 90fm/c, respectively, at the energies of Ecm = 0.1, 4 
and 25 MeV. The backward evolution initially follows very well the forward evolution with 
suppression, especially well at the two higher energies. Gradually, however, deviations build 
up, and the system extents evolved in the two ways, with and without suppression, separate 
at all energies and for every case for times prior to the initial compression. This is indicative 
of the fact that the system, at either of the studied energies, does not evolve back into 
the initial state. A particularly dramatic situation develops at the lowest of the energies, 
where the system evolved backward in time stays compact and oscillates in size, rather than 
separating into the two original fragments, for which system extent would increase with 
decreasing time. At the two other energies, we can observe a reduced slope in the dependence 
of system extent around the initial time. That indicates a reduced speed and/or a reduced 
mass for the fragments, compared to the original initial state, both implying a certain level 
of internal excitation. Clearly, with element suppression, the time reversal symmetry is 
broken at the level of global system characteristics. The degree of violation of time reversal 
symmetry, due to element suppression at fixed xq, depends on the collision energy. This is 
similar to the situation with the impact of element suppression on the forward time evolution. 
Further insight into the time-reversal violation induced by off-diagonal element suppres- 



sion is provided by Fig. 27, that compares the initial system density to the density at t = 
obtained from evolving the system first forward and then backward in time, at the three 
collision energies. For the system evolved forward and backward at the lowest energy of 



Ecm = 0.1 MeV, with element suppression, we can directly see in Fig. [27] a fused system at 
t = 0, that has remained fused and oscillating since the initial contact. At Ecm = 4 MeV, 
with element suppression, the system evolves back into two highly excited fragments around 
t = 0, that are slowed down compared to the original initial state. The level of excitation 
is evidenced by the depth of the dips in density for these fragments in the central panel of 
Fig. |27| Finally, at Ecm = 25 MeV the system evolves back into two excited fragments with 
a low-density neck formed in-between. 

At each energy, the system that is evolved back to t = 0, with suppressed elements in the 
density matrix, acquires some features typical for the final states at the specific energy, be 
that fusion, fragment excitation, or formation of a neck. This indicates that the correlation 
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patches within the density matrix that may be far from the x = x' axis late in the system 
forward evolution, and that move back towards the axis when the time is reversed, are 
essential in the restoration of the density matrix of the coherent initial state. When these 
off-diagonal elements are erased, the initial state coherence is lost and it is not possible to 
recover a well separated, two-slab system. Note, also, that for the forward time direction, 
the entanglements are of far less importance because coherence is less of an issue at late 
stages of nuclear reactions. 

One aspect of the collisions that has not been addressed yet is their momentum content. 
In the next Section, we remedy this by analyzing a mixed, configuration and momentum 
space, representation of the density matrix. The discussion will also bring in another inter- 
pretation for the erasing of the far off-diagonal elements within the spatial representation. 



IX. WIGNER DISTRIBUTION 



As has been discussed in the Introduction, high-energy central reactions are commonly 
described in the literature in terms of the semiclassical BE. The dynamic quantities in BE- 
type of approaches are the so-called Wigner functions, distributions of particles jointly in the 
configuration and momentum space. Wigner distributions have been also employed in the 
literature in the context of TDHF, in particular to gain guidance in supplementing TDHF 
with a collision term [3114371 139] . Within nuclear structure, the Wigner functions represented 
a departure point for semiclassical approximations [2Qj. Elsewhere, distributions in the 
mixed configuration and momentum representations have been found to provide interesting 
insights into aspects of quantal dynamics [55]. For ID systems, as will be evident, the 
Wigner functions are particularly straightforward to visualize. 

The Wigner function jy^ results from Fourier transforming the density matrix in the 
relative coordinate: 

fw{x,p) = dxre~^p{x + Xr/2,x — Xr/2). (46) 

The transformation yields a real function fw that, in the classical limit, can be interpreted 
as a phase-space distribution [55]. This function is, however, generally not positive definite. 
In quantum mechanics, it only acquires that property after fw is averaged by folding it with 
a sufficiently wide Gaussian in position and momentum. Consistently with its interpretation. 
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projecting the Wigner function onto the x-axis yields, up to a factor, the ID spatial density: 

dp 1 

fwix,p) = p{x,x) = -ni{x) . (47) 

Similarly, projecting onto the p-axis yields, up to a factor, the density in momentum space. 

Under proper conditions [2], which include meeting in different ways the semiclassical 
limit, the KB equations yield a BE for /^y. The collision term in BE is associated with the 
correlation contributions to the KB equations, cf. the l.h.s. of Eqs. ([s]) and Q. This collision 
term, in particular, includes the effects of statistics. Aside from being highly exploited in 
central high-energy nuclear reactions t26J, BE with this type of a collision term has been 
beneficial in other fields O [50] . 

When the semiclassical limit is applied on top of the mean-field dynamics, the Vlasov 



equation for fy^ emerges from the equation of motion (10) for the density matrix. This 



represents a reduced form of BE with effects of collisions suppressed. Specifically, given the 
peaking of the density matrix in the relative coordinate around zero, we can attempt to 



expand the difference of mean- field potentials on the r.h.s. of Eq. (10): 

dC/ 



dx 



(48) 



Combining this approximation with a Fourier transformation in the relative coordinate of 



both sides of Eq. (10), yields the Vlasov equation for fyi/: 

dfw , P dfw dU dfv 



+ 



'w 



0. 



(49) 



dt m dx dx dx 
The Vlasov aspect of the equation is in the self-consistent evaluation of U from the single- 
particle density. The equation for the evolution of fw under the influence of an external 
potential U is of an identical form. 



Following Eqs. (32) and (33), the width of the main peak around the diagonal of the 



density matrix for our ground-state slabs is about ttH/pfi ~ 3.2 fm, see also Fig. 10 For 



a general potential [/, the validity of Eq. (49) requires a slow variation of the potential over 



that specific distance. As the momentum scale is compared with the scale in position and h 
is involved, this condition is semiclassical in nature. Notably, for a U that is largely linear 



in density at low to moderate densities, cf. Eq. (17), the condition of slow variation will not 



be satisfied for our slabs, given their shape, see Fig. 2 On the other hand, Eq. (48) may 



yield coarsely correct results. Otherwise, since the approximation of a linear expansion of U 
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at best addresses the immediate vicinity of the diagonal in the density matrix, the issues of 
any structures moving away from the diagonal are left out in the Vlasov description. Note, 



however, that in the case where U is quadratic in position, Eq. (48) is exact. For such 
a potential, it does not matter whether there are any off-diagonal structures present or not. 
In particular, the Wigner function for the HO initial state of the adiabatic evolution will 



represent a static solution to Eq. (49). 



In spite of these limitations, the Vlasov equation (49) is still useful as a reference when 
analyzing collisions in terms of Wigner functions. The equation can be written in different 
ways and, in particular, can be cast into the form of a single-particle Liouville equation. 
Upon defining the single-particle energy as 

2m 



e{x,p;t) = ^ + U{x,t), (50) 



the Vlasov equation may be rewritten as 

dfw , de dfw de dfw dfw 9 f de \ d f de 

+ 7, 7^ 7, = -ITT + ^\ ^JW]+Tr \ -TT JW 



dt dp dx dx dp dt dx \dp J dp \ dx 

_ dfw 



dt 



+ tCfw = 0. (51) 



A combination of the terms on the l.h.s. of Eq. (51) may be recognized as the Poisson 
bracket of the Wigner function with the single-particle energy. Upon rewriting of the l.h.s. 



of the Vlasov equation after the first equality in (51), the Vlasov equation acquires the form 



of a continuity equation in phase space, with — ^ j representing the local phase-space 
velocity- vector and (^^fw,~^fw^ representing the phase-space flux. The semiclassical 
limit of the Liouville operator in phase space is then given by 

d fde \ d fde \ 



dx \dp J dp \dx 
If the equations for the phase-trajectories {x{t),p{t)) are introduced, 

dx de , dp de 

Tt = dp d^ = -a^' 

the Liouville equation is found to imply that the phase-space density is constant along the 
trajectories, 

j^fw{x{t),p{t)-t) = 0, i.e. fw{x{t),p{t)-t) = /h.(x(0),p(0);0). (54) 
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The classical phase space density described by the Vlasov equation is just the number of 
particles per phase-space volume. Within the mean-field dynamics, those particles move 



in time along the trajectories given by (53). The conservation of the density along these 
trajectories may be associated with the fact that the classical evolution, represented by 



Eqs. (53), preserves the phase-space volume spanned by the trajectories. The latter follows 
from the fact that the phase-space divergence of the phase-space velocity vanishes and it 
may be also determined by examining directly the changes in the phase-space Jacobian: 

if'fi-f>^0. hence f'fj-f '=1. (56) 

dt d{x{0),p{0)) ' d{x{0),p{0)) ^ ' 



The first equality follows from employing the equations of motion, Eq. (53 ), in calculating the 
derivative of the Jacobian. The practical outcome is that both the phase-space distribution 
and the phase-space itself behave as an incompressible fiuid. In a sense, the distribution 
may be viewed as a set of markings on locations within a phase-space fiuid that evolves with 



time in an incompressible manner, consistent with Eq. (53). 

In the following, we shall examine the Wigner function for the collision of A = 8 slabs at 
Eqm = 25 MeV. We choose that particular energy because the Wigner function is then most 
spread out in momentum space, making many details particularly visible. With calculations 
carried out in a periodic box of —L < x < L, our Wigner function is defined for momenta 
mTh/2L, where n is integer. With L = 25 fm, this provides a momentum resolution of Ap = 
7rh/2L = 12.4MeV/c or, equivalently, a wavevector resolution of Ak = n/L = 0.063 fm~^. 



Note that the integration over Xr in (46 ), for the periodic box, extends over a distance of 2L. 



The left panel in Fig. 28 shows the Wigner distribution for the initial state of our system 
at Eqm = 25 MeV. With the initial density matrix equal to the sum of density matrices for 



the two slabs, cf. Eq. (30), the Wigner function is a sum of the Wigner functions for the 



two slabs, each represented by a ring structure in Fig. [28j The centers of those rings are 
displaced from the origin by the displacement of the slabs from x = (their centroids lie at 
±x = 7.5 fm, cf. upper left panel of Fig. [7| and by ±P/h = 1.10 fm~^ from p = 0. The latter 
displacement follows from the fact that multiplication of a density matrix by a phase-factor 



Q±tPxr^ cf. Eq. (29), leads, under Fourier transformation, to a displacement of the Wigner 
function in momentum by ±P. Relative to its center, a ring represents the Wigner function 
of a ground-state slab at rest. 

The Wigner function of each slab is basically a sum of the Wigner functions for each filled 
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orbital, cf. Eqs. (13) and (46). For orbitals that are parity eigenstates, the orbital Wigner 
functions are invariant under inversion in (x,p)-space. Given that the density of the ground- 
state slabs is similar to that in the ground state of the corresponding HO potential, cf. Fig. |2| 
we expect similarity in the Wigner functions as well. Classical HO orbits for constant single- 
particle energy, e = p^/2m + are ellipses in phase space. They would become 
circles if the momentum axis was rescaled by a factor of 1/ (mfi) and, consistently with 
this, the Wigner functions of HO orbitals become isotropic with such a rescaling. They 
are given by Laguerre polynomials multiplying Gaussians in distance from the origin 
As the lowest orbital is a Gaussian, the corresponding HO Wigner function is also a Gaussian, 



cf. Eq. (46). Clearly, however, a Gaussian does not explain the ring structures in the left 



panel of Fig. 28 The responsible party for such a structure is actually the second orbital, 
which is antisymmetric and, thus, has a node at the origin. For the HO potential, the ring 
may also be imagined in terms of a collection of classical elliptical orbits. With time, the 
particles circle on those orbits in the clockwise direction with a period of 2Ti/Vt ~ 92fm/c, 



cf. (B4). For a mean-field potential, the bounded orbits will not be substantially different 
and therefore they show an elliptical structure. 



With the single-particle energies of Eq. (50) being continuous functions of location in 



phase space, an evolution consistent with the classical limit of Eq. (49) should preserve the 
topology of the original Wigner function. Consequently, if the original topology is in the form 
of two rings, those two rings should be recognizable later in the evolution. Large distortions 
in the shape might appear, but their topology should be preserved in the time evolution. 
Consequently, the same set of values for the function inside and along the circumference of 
the rings might be expected. 



The next stage of the reaction, represented in the central panel of Fig. [28} presents the 
actual quantum evolution of the Wigner function and corresponds to the stage of overlap 
and maximal compression of the two slabs. In spite of the maximal compression, which 
corresponds to a single compound slab in configuration space, the two rings in the Wigner 
distribution can be perfectly distinguished. In the initial state, the two rings are separated 
in phase space and this topology is preserved here, even though the rings overlap in con- 
figuration space. The rings are, however, distorted compared to the initial state. After the 
slabs come into contact, the nucleons that get first affected are those that move towards the 
other slab. With the potential barrier between the slabs gradually disappearing, the fastest 
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nucleons will move into the other slab. Slower nucleons, being in the course of getting re- 
flected from the potential barrier of their parent slab, will end up losing their momentum 
only partially at the merger stage of the slabs, as the barrier disappears. Finally, the slowest 
nucleons, moving towards the outer boundary of their parent slab, continue to do so at that 
stage. This produces a characteristic trapezoidal shape for each of the rings in the center 



panel of Fig. 28, representing t = 30fm/c. 

In spite of the similarities between the initial state and the maximal compression Wigner 
functions, there are structures in the latter functions that are of a pure quantum nature. One 
nonclassical aspect of the phase-space distribution is that of the thinning of each ring around 
a line starting from the side of the ring lowest in momentum magnitude and ending at the 
ring's center. The values of the distribution drop there everywhere radially across the ring, in 
contradiction to what is allowed for a classical distribution starting from that in the left panel 



of Fig. 28 Another nonclassical aspect are the values of the distribution right in-between 
the two rings. For a classical evolution, we would expect a valley with vanishing values 
of the Wigner function. Instead significant negative values are found there. The specific 
valley has been observed before in ID slab collisions [36l |39] and it has been discussed in 
the context of the lack of collisions in TDHF. Correlations, if included, would generally act 
to thermalize the single-particle Wigner function, gradually filling the valley j37j . 

The Wigner distribution for the late stage of the reaction, t = 80fm/c, is shown in the 



third panel of Fig. 28 The distribution is complicated and we may have a hard time recog- 
nizing traces of the original rings. In large part this is due to quantum-mechanical effects 
whose presence could be anticipated given the substantial off-diagonal structures encoun- 
tered in the spatial representation of the density matrix. In spite of those complications, 
it becomes generally apparent, from the sequence of phase-space images in Fig. [28} that 
the two fragments leading in the configuration space, that emerge from the collision, see 
also Fig. [7} represent the fastest nucleons from the opposing slabs. Those nucleons have 
punched through the matter and have remained the fastest within the system as the col- 
lision progressed. The original ring structures may be still recognized at t = 80fm/c in 
the regions immediately adjacent to the peaks in phase space that represent the leading 
fragments. From there on, towards lower momenta, streaks may be observed that can be 
attributed to the same original slab as the leading fragment. However, the ring structure 
cannot be recognized within the respective streak region. In fact, three or more parallel 
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ridges of greater intensity may be counted there within the remnants of each of the slabs, 
while the classical dynamics could explain the presence of only two parallel ridges of spe- 
cific spatial length, formed over a specific time. A closer examination reveals many more 
fainter and more narrowly spaced parallel ridges aX t = 80fm/c, down to the momentum 
resolution from the inverse size of the calculation region. These clearly represent interfer- 
ence patterns that mar any phase-space structures that could be attributed to a classical 
dynamics. Between positive values of the function for the ridges, negative values occur in 
the valleys. The gross pattern of the distribution at t = 80fm/c is that of a developing 



Hubble correlation between position and average momentum, cf. Eq. (38). 

In the preceding Sections, we have discussed the role of far off-diagonal elements of the 
density matrix and the effects of a potential suppression of those elements. We will now 
discuss the interrelation between the suppression of elements and the momentum space 



features of the Wigner function. Clearly, given Eq. (46), the Wigner function must be 
affected by element suppression. To gain an understanding of the effect, let us denote as 
fw the Wigner function for some density matrix, p. Let us further denote as fp the Wigner 
function for a matrix obtained from p by suppressing the far-away elements according to 
some profile function, V{xr): 

r ipxr 

fp{x,p) = / dxr V{xr) p{x + Xr/2, X — Xr/2) , (56) 

where V{0) = 1. In the case that we have discussed in Sec. VI, P is generated by the 
exponential suppression factor associated to the superoperator field, W. For a sharp cutoff, 
with matrix elements put to zero at \xr\ > Xq, the profile function would be V{xr) = 0{xq — 
\xr\)- As generally known, a product in configuration space transforms into a convolution 
in the conjugate momentum space. Upon expressing the density matrix p in terms of fw, 



following Eq. (46), we find: 

fv{x,p) = J dp'Vip - p') fwix,p') , (57) 

where 

^(P) = ^ I e-"^ V{Xr) . (58) 
Given that V{xt = 0) = 1, we find that 

dpV{p) = l, (59) 
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and we therefore observe in Eq. (57) that the Wigner function /-p amounts to an average 



out in momentum of the function fw The weight function in averaging is the Fourier 
transform of the profile function. In the case of the exemplary sharp cutoff profile above, the 
weight function is, after (58), V{p) = :^sm^, with the range of averaging in momentum 
of ~ TTh/{2xo). The general conclusion that follows from this investigation is that the 



suppression of far-away elements should lead to a blurring of finer momentum-structures 
within the Wigner distribution. 



The above expectation is tested in Fig. 29, that shows the Wigner distribution at the 
late stage, t = 80fm/c, of the Ecm/^ = 25MeV collision, from calculations with different 



cutoff parameter xq in the superpotential W of Eq. (43). At this stage of the collision 



the Wigner function is particularly nuanced. If we compare the left panel in Fig. 29 to the 



right panel in Fig. [28} representing the standard evolution, we can see that even the cutoff 
at Xq = 15 fm significantly blurs the Wigner function. As the cut-off in W is reduced, from 



the left towards the right panel in Fig. 29, the Wigner function gets progressively more 



blurred, with the finer interference patterns disappearing. The patterns that remain may 
be actually attributed to the abruptness in suppressing elements as a function of Xr- For 



comparison, we further show in Fig. 30 the Wigner function from the standard evolution, 



i.e. that represented in Fig. 28, but now with the function directly blurred at t = 80fm/c by 
folding it in momentum with a Gaussian profile function of a prescribed width. Specifically, 
we make the momentum range for the profile function. 



V{p) 



a 



exp 



'2^h ' V 2/12 

about the same as that for the weight from a sharp cutoff at xq, with 



(60) 



a 



2x0 

TT 



(61) 



In comparing Figs. 29 and 30, we see a good degree of agreement between the corresponding 
blurred distributions, for the two larger values of Xq, in the regions of significant values of 
the distribution. For xq = 5fm, the agreement is qualitative. Note that the convolution in 



momentum with Eq. (60) does not change the density for the standard evolution and the 



densities for evolutions with and without cutoffs have been already compared in Fig. 23 



Several conclusions, pertaining to the single-particle quantities and their evolution, may 
be drawn from the above considerations. Suppressing the far-off diagonal elements in the 
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density matrix is equivalent to smoothing the Wigner function in momentum. The more se- 
vere the suppression of elements, the greater is the level of smoothing of the Wigner function. 
This equivalence is important for possible extensions, where correlations are incorporated 
into the single-particle approach using the KB equations. In that case, NGFs with different 
time-arguments need to be considered, see Sec. [TT| The smoothing of Wigner distributions 
for the purpose of discarding redundant information might be generalized in a more straight- 
forward manner in the correlated case than the suppression of off-diagonal matrix elements. 
We can further observe that, for moderate levels of smoothing, it is sufficient to follow only 
the self-consistent evolution of the smoothed Wigner distribution, to arrive at a faithful 
information on the smoothed distribution in the future. Finally, it is important to note 
that, in the region of significant values of the Wigner distribution, the particular details of 
smoothing, other than its overall range, are likely to be irrelevant, as seen by comparing 



Figs. 29 and 30 



X. CONCLUSIONS 



In spite of their considerable general potential in describing quantal time-dependent sys- 
tems, NGF techniques have been, so far, underutilized in nuclear physics. Specific potential 
applications within nuclear physics include the description of large-amplitude collective mo- 
tion and of central nuclear collisions. The techniques could yield a seamless description for 
these areas and for nuclear structure [19]. The scarce use of these techniques to study nu- 
clear reactions in practice could be attributed to the anticipated large computational effort. 
In this paper, we have investigated whether all the information that would be naturally 
followed within the NGF approach, actually needs to be maintained in describing nuclear 
reactions. This is the major source of computational effort within the NGF method and 
suitable information suppression techniques would therefore have the potential to make such 
calculations feasible. 

For this purpose, we have studied ID slab collisions within the mean-field approxima- 
tion. While our starting point has been equivalent to the TDHF approach, the manner of 
discarding of information pursued here could not have been investigated within an approach 
relying on wavef unctions only. We have derived correspondences between 3D and ID systems 
that help to comprehend the features of the latter. In solving the single-particle evolution. 
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we have demonstrated in practice the possibihty of arriving at a mean-field ground-state, 
through an adiabatic transformation of the interaction. We hope to be able to arrive sim- 
ilarly at the ground state for correlated dynamics. Following our general goal, within the 
mean-field approximation, we have examined the structure of single-particle density ma- 
trices for collisions. The matrices are generally peaked along the diagonal in their spatial 
arguments. For the ground state, the width of the peak of significant values of the matrix 
elements around the diagonal is proportional to the inverse of the Fermi wavevector. Other- 
wise, initial density matrices are confined within square-like structures in the {x, x') plane, 
with the side of their compact support equal to the spatial extent of the respective slab. 

In the dynamical evolution associated with a reaction, significant values of far ofi^-diagonal 
elements may appear due to the entanglement of single-particle wavefunctions for nucleons 
that have the possibility of moving out with different momenta and different fragments from 
the reaction. Applying an absorptive superpotential within an evolution of the density ma- 
trix, we have demonstrated that the entanglement patches in the far-off diagonal regions of 
the matrix have little effect on the evolution close to the diagonal. This can be attributed 
to the fact that the emerging fragments have little chance to encounter each other again in 
an expanding system, particularly at low relative velocities where phases in the wavefunc- 
tions, recorded in the patches, might matter. Within the structure of the density matrix, 
those entanglement patches generally travel away from the diagonal of the matrix, never to 
come back. Even when the element suppression with the superpotential is severe and the 
density matrices of the individual fragments are modified, we find that the intrinsic collective 
motions for the fragments are only affected at a semiquantitative level. Also, we find that 
the energy conservation is quite robust when off-diagonal elements of the density matrix 
are suppressed. These findings bode well for the possibility of carrying out realistic NGF 
calculations of nuclear collisions, where the discarding of information, such as in elements 
far away from the diagonal of function arguments, becomes an absolute necessity. On the 
other hand, we find that maintaining the far off-diagonal elements is important for pre- 
serving time-reversibility of the mean-field evolution. Under time reversal, the correlation 
patches within the density matrix travel back towards the diagonal of the matrix, where 
they contribute to the restoration of original density matrix for the initial state. If the cor- 
relation patches are removed, the state restored in the course of backward evolution acquires 
features that are actually characteristic for final states of the collisions within the specific 
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energy regime. 

As is well known, Fourier transformation of the density matrix in relative coordinates 
yields the Wigner function. The semiclassical Vlasov equation for the Wigner function re- 
sults under the assumption that the density matrix is peaked sharply enough around the 
diagonal within its arguments so that the mean-field in the equation of motion can be ex- 
panded in the distance away from the diagonal. In aiming at suppressing the off-diagonal 
elements of the density matrix, we do not insist on the validity of the semiclassical ex- 
pansion. Rather, we aim at a situation where the element suppression is dialed to achieve 
a compromise between the computational effort to be undertaken and the detail needed in 
the description of physical processes. The proximity of a system to the semiclassical limit 
would then just act to improve the chances for reaching a satisfactory compromise. The sup- 
pression of the far-away elements for the matrix, in the spatial representation, is equivalent 
to smoothing of the matrix within the momentum variable in the Wigner representation. 
Such a smoothing procedure removes fine fringe patterns in the structure of the Wigner 
function. The equivalence between off-diagonal element suppression and momentum space 
smoothing is also important for the way in which discarding of information needs to be 
generalized to the case of NGF with correlations. In that case, the smoothing might be 
a more natural generalization when NGF with different time arguments are considered. 

To our knowledge, this has been the first time when the off-diagonal elements in the 
density matrix have been investigated within the context of nuclear collisions. The original 
scope of this investigation has been modest, aiming at a reduction of the computational effort 
within the NGF approach, with the testing ground being the ID mean-field dynamics. In the 
end, however, the investigation touches upon some general aspects of quantal many-body 
dynamics. Possibly the most important of those is the growing redundancy of information 
within a system undergoing expansion. 
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Appendix A: 3D interpretation of ID densities 

Results of ID calculations may be interpreted in terms of 3D matter, when assuming 
that the 3D matter is uniform in two directions y and while generally nonuniform in 
the remaining x direction. In the two uniform directions, the matter is described in terms 
of a set of plane- wave wavefunctions that multiply the wavefunctions labeled with a, de- 
scribing variation in the x direction and incorporating spin and isospin degrees of freedom. 
The transverse set is assumed frozen, independent of a or time. Under these assumptions, 
the density of 3D matter, n^{x,t) = n{x,t), becomes proportional to the density, ni{x,t), 
of ID matter calculated from the ID density matrix or its wavefunctions (l)a{x,t), 

n{x,t) ^ ^ni{x,t) . (Al) 

The scaling factor ^ may be established following the requirement that the energy of 
uniform matter needs to minimize at the normal density no- In the derivation, we shall 
invoke the Hugenholtz-van Hove theorem, both for the standard isotropic 3D matter and 
for the ID matter to be interpreted in the 3D manner. The Hugenholtz-van Hove theorem 
states that the energy per nucleon and the Fermi energy should coincide when the system is 
in its ground state. In the mean-field approximation, both the energy per nucleon and the 
Fermi energy consist of kinetic and mean-field terms. The Fermi energy at normal density 
is 

SF1,3 = ^ + ^(^o) , (A2) 

where the Fermi momentum, pp, is either that of ID matter or that of the standard 3D 
matter. The total energy per nucleon is 

for the 3D matter and 

for the ID matter. Here, eu is the contribution to the energy associated with the mean-field. 
In the ID case we consistently include only the kinetic-energy contribution associated with 
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the X direction. The Hugenholtz-van Hove theorem states that 



£^F3 = es . 



(A5) 
(A6) 



Upon subtracting the equations side-by-side, the mean-field contributions drop out and we 
find the relation 



Pfi 



at normal density. We further have 



rii 
no -- 



PF3 



^Pf3 



(AT) 



(A8) 
(A9) 



Upon combining Eqs. (Al), (A7), (A8) and (A9), we find 



vrn, 



2\ 1/3 




3 V 6z/2 



1.04 



no \ 2/3 
u J 



(AlO) 



Thanks to Eq. ( Al ) and the considerations that follow ( Al ), 3D mean-field parametrizations 



from the literature can be immediately adapted for ID. 

The considerations above can further provide an estimate for the thickness of a slab for 
a given A in the ID interpretation: 



A 

ni 



no 
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3 I, 6z/2 no 



1/3 



A ~ 0.80 fm A. 



(All) 



In the above, we have employed no ~ 0.16 fm~^. Otherwise, when we quote a mass number 
A in the ID interpretation, this represents the number of nucleons in the 3D interpretation 
within a transverse area of 1/.^. 



Appendix B: Oscillator frequency for ID slabs 

To minimize the adjustments that the system needs to go through during the adiabatic 
switching on of the interactions, one needs to choose suitable precursor HO states. These 
states are characterized by two parameters: the total number of filled shells, Ns, and the 
oscillator frequency, Q. The first parameter determines the number of nucleons A in the 
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ID interpretation of the results and, further, the number of nucleons per unit area in the 



3D interpretation, see Sec.|ITT]and Appendix |AJ For spin-isospin saturated systems, studied 
here, A = u Ng, with z/ = 4. To achieve a density in slab interior which is close to the normal 
density for uniform matter, the frequency Q needs to be adequately correlated with A. 

To establish what Q should be used for a given A, we shall follow a procedure analogous 
to that employed in the literature for the 3D HO. The well-known 3D result of HQ = 
4:1A~^^^ MeV [29] is found by combining the virial theorem for HO with an expression for 
the mean square radius of a uniform sphere. Correspondingly, we start out from the virial 
theorem for a single-particle state a of the ID HO: 

-mQ' {x')^ = — + , (Bl) 

where {x'^)^ is the mean square position within the state a and is the oscillator 
quantum number. If A^^ states are filled, the corresponding quantum number values are 
no, = 0, 1, ■ ■ ■ ,Ns — 1. The mean square position within a slab is 



Upon summing up both sides of Eq. (Bl) over the A^^ shells, we find 

- ^''^^^ - ^'-^ rR3i 

2m{x^)A~ 2mu{x^) ' ^ ' 

where we have employed A = vN^. For a uniform slab of total thickness £, we would have 
had (x^) = £^/12. If we aim, for large A, for a slab of average density A/^ = rii, we need to 
choose an HO frequency of 

_ rij _ nl 108 MeV 
m uA m pA A ' 

where we have made used of the connection between ID and 3D density derived in Ap- 
pendix |X} The utihty of the derived relation, for adiabatic switching, can be tested a pos- 
teriori by comparing the density profiles of the predecessor HO state and that of the cor- 
responding self-consistent slab. In the example provided in Fig. |2j the profiles of both the 
precursor HO slab and the self-consistent mean- field one are very close. We have observed 



a similar agreement over a wide range of particle numbers A, which suggests that Eq. (B4) 
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indeed yields a good choice of the frequency for the starting state. 
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FIG. 1: (Color online) Time evolution of the energy per particle (upper panel) and the size of the 
slab (lower panel) when starting from an HO configuration with Ns = 2 filled shells (or A = 8 in 
ID interpretation) at time t = —1000 fm/c and transforming the single-particle potential to the 



mean-field form, according to Eq. (23) with Eq. (25) and (26). Different values of the transition 
time T are considered. Inset into the top panel shows a magnified portion of the time evolution 
of the energy. For reference, the mean-field results from static Hartree-Fock solution are shown as 
straight solid lines. 
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Adiabatic switching: local density 




I 0.05 



FIG. 2: (Color online) Evolution of the density profile, when starting from an HO configuration 
with Ng = 2 filled shells (or A = 8 in ID interpretation) at time t = — lOOOfm/c and transforming 



the single-particle potential to the mean- field form, according to Eq. (23) with Eq. (25) and (26), 
r = 40fm/c and tq = — 600fm/c. 
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FIG. 3: (Color online) Time evolution of the energy per particle (upper panel) and the slab 
width (lower panel), when starting from an HO configuration with Ng = 2 filled shells at time 
t = —1000 fm/c and transforming the single-particle potential from the HO to the mean-field form, 



following different switching functions: of the Fermi-Dirac form [see Eq. (26)], linear and piecewise 
quadratic. 
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FIG. 4: (Color online) Evolution of the CM density profile, represented by solid line, for a collision 
of two A = 8 slabs at the collision energy of Eqai/A = 0.1 MeV. The dashed line in the last panel 
represents the density profile for the ground state of an ^ = 16 slab, provided for comparison. 
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FIG. 5: (Color online) Time evolution of system size for a collision of two ^ = 8 slabs for different 
indicated collision energies. Oscillating sizes in the late stages of the evolution represent fused 
systems. Sizes that grow monotonically correspond to systems that did not ultimately fuse. 
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FIG. 6: Evolution of the CM density profile for a collision of two A = 8 slabs at the collision energy 
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FIG. 8: (Color online) Schematic illustration showing the origin of the characteristic features in 
the off-diagonal structure of a single-particle density matrix. In the case of localized single-particle 
states (left panel), the density matrix is confined to a square-like region. Fragmented single-particle 
states (right panel) give rise to a patch structure for the density matrix, with patches extending 
arbitrarily far away from the diagonal. 
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FIG. 9: (Color online) Intensity plots representing the real (upper panels) and the imaginary part 
(lower panels) of the scaled density matrix, p{x , x' ; t) , for a collision at Ecm/A = 0.1 MeV. 
The left, central and right panels show the initial, overlap and late stages of the collision, respec- 
tively. The scaling aims at making the values along the diagonal coincide with the 3D density 
shown in Fig. |4j 
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FIG. 11: (Color online) Same as Fig. [o] but for a collision at Ecm/A = 4 MeV. 
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FIG. 12: (Color online) Same as Fig. 10 but a the collision at Ecm/A = 4 MeV. 
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FIG. 14: (Color online) Same as Fig. 10 but for a collision at Ecm/A = 25 MeV. 



68 



X' 




FIG. 15: (Color online) Structure of the superpotential W, and of the corresponding suppression 
of matrix elements, within the {x,x') plane. In a band around the diagonal, \x — x'\ < xq, 
darkened in the figure, the superpotential vanishes and the elements are not suppressed. At 
separations xq + do < \x — x'\ < L, shown as white, the superpotential is strong and the matrix 
elements are strongly suppressed. The superpotential and suppression change gradually in-between. 
The periodicity of the system gives rise to repeated bands of directly unaffected and suppressed 
elements in the corners of the computational space. 
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FIG. 16: (Color online) Time evolution of the energy per particle (upper panel) and of the system 
extent (lower panel) in a collision at Ecm/^ = 0.1 MeV for different values of the cutoff parameter 
xq associated to the suppression of off-diagonal elements in the density matrix. 
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FIG. 17: (Color online) Density profiles in the Ecm/^ = 0.1 MeV collision, at sample times, 
for different values of the cutoff parameter xq. The profiles are staggered by 0.1 fm^'^ to provide 
better insight into details. From bottom to top, the lines are for the standard evolution, and then 
Xq = 20, 15 and 10 fm, respectively. 



70 



^ 0.4 
0.3 h- 

>< 

0.2 
X 0-1 



q 



t=0 fm/c, X=7.5 fm t=300 fm/c, X=0 fm 
T — I — I — I — I — I — I — I — r 




J I I I I I I I L 



1 1 


1 1 1 1 1 1 
i '! 

!'.'• - 






— 


.'^ r ~ 


— " ^ 

\ 

\ 


A A/- 






■/•• 103 



t=500 fm/c, X=0 fm 
n — I — I — I — I — I — I — I — r 



/ 1 
I i 



I 1 
1 1\ i 




J I I I I I I I L 



0.4 ^ 



-20 -10 10 20 -20 -10 10 20 -20 -10 10 20 



X 

>< 

0.2 

>< 

0-1 >< 



q 
> 



X [fm] 



[fm] 



X [fm] 



FIG. 18: (Color online) Real part of the scaled density matrix p{x,x' ]t) for the slab collision 



at Ecm/A = 0.1 MeV, along a line x + x' = 2X = const, perpendicular to the diagonal x 



X 



in the {x,x') plane, for different values of the cutoff parameter xq. The results are shown against 
the coordinate difference x — x' = Xr, &t sample values of X and t. The results for different xo 
are staggered by 0.1 fm^^ to provide better insight into details. From bottom to top, the results 
are for the standard evolution, and then xq = 20, 15 and 10 fm, respectively. The locations where 
element suppression sets in, \xr\ > xq, are marked with triangles. 
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FIG. 19: (Color online) The same as Fig. 16 but for a collision at Eqm/-^ = 4 MeV 
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FIG. 20: (Color online) The same as Fig. 17 but for a collision at Eqm/A = 4 MeV 
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FIG. 21: (Color online) The same as Fig. 18 but for a collision at Ecm/A = 4MeV 
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FIG. 22: (Color online) The same as Fig. 16 but for a collision at Ecm/A = 25 MeV. 
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FIG. 23: (Color online) The same as Fig. 17 but for a collision at Eqm/-^ = 25MeV. Note the 



difference in vertical scales for the left panel and the two remaining panels. 
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FIG. 24: (Color online) Intensity plots representing the real part of the scaled density matrix, 
(,1^ p{x , x' ; t) , for a collision at Ecm/A = 25MeV at t = 80 fm/c in calculations with different off- 
diagonal cuts, xq. The corresponding intensity plot for a calculation without element suppression 
may be found in Fig. Il3l 
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FIG. 25: (Color online) The same as Fig. [21] but for a collision at Ecm/A = 25 MeV. Note the 
differences in the vertical and horizontal scales for the left panel and the two remaining panels. 
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FIG. 26: (Color online) Evolution of the system extent for the collision of two ^ = 8 slabs at 
three energies. The standard, fully reversible evolution is represented by dotted lines. The forward 
evolution, with the xq = 10 fm suppression, is represented by solid lines. The backward evolution, 
with the suppression, is represented by dashed lines. In each case, the reversal of evolution is 
applied at the end of the displayed time interval. 
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FIG. 27: (Color online) Comparison of the initial density (solid lines) to the density obtained at 
t = when first evolving the system forward and then backward in time, with element suppression 
in the density matrix, for a collision of two A = 8 slabs at three energies. 



t=0 fm/c t=30 fm/c t=80 fm/c 




X [fm] X [fm] X [fm] 



FIG. 28: (Color online) Intensity plot of the Wigner distribution, fw{x,p), for the collision of 
^ = 8 slabs at Ecm/A = 25 MeV. The left, center and right panels represent, respectively, the 
sample times of t = 0, 30 and 80 fm/c. The vertical scale is A: = p/h. 
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FIG. 29: (Color online) Intensity plot of the Wigner distribution, f\y, att = 80fm/c, for Ecm/A = 
25 MeV A = 8 slab collisions, from calculations with different element suppressions. The left, center 
and right panels represent, respectively, calculations with the cut-off parameter value of xq = 15, 
10 and 5fm. The vertical scale is fc = p/h. 
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FIG. 30: (Color online) Intensity plot of the Wigner function, /p, obtained by convoluting the 
function fw from the standard calculation of the Ecm/-^ = 25 MeV collision, for t = 80fm/c, with 



a Gaussian [see Eq. (60)]. The width of the Gaussian for the three displayed panels has been taken 



in proportion to the cutoff, xq, in the three respective panels of Fig. 29, cf. Eq. (61). With this, 
there is a match between the range of Gaussian averaging and the expected range of blurring of 
the Wigner function due to the suppression of elements in the evolution of the density matrix. 
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